^D-fli34  830 THE  DEVELOPHENT  OF  fl  MUMERICflL  SOLUTION  TO  THE 

TRANSPORT  EQUATION  REPORT  2  COMPUTATIONAL  PROCEDURES 
(U)  COASTAL  ENGINEERING  RESEARCH  CENTER  VICKSBURG  MS 
UNCLASSIFIED  R  A  SCHMAL2  SEP  83  CERC-83-2-2  F/G  8/3 


VF3a 


MISCELLANEOUS  PAPER  CfcRC-83-2 


THE  DEVELOPMENT 
OF  A  NUMERICAL  SOLUTION 
TO  THE  TRANSPORT  EQUATION 

Report  2 

COMPUTATIONAL  PROCEDURES 


by 

R.  A.  Schmalz,  Jr. 

Coastal  Engineering  Research  Center 
U.  S.  Army  Engineer  Waterways  Experiment  Station 
P.  O.  Box  631,  Vicksburg,  Miss.  39180 


September  1983 
Report  2  of  a  Series 

Approved  For  Public  Release;  Distribution  Unlimited 


>- 

Q- 

O 


PTIC 


'V;.  NOV  1  3  1983  ;;j 


Prepared  lor  U.  S.  Army  Engineer  District,  Mobile 
Mobile,  Ala.  36628 

11  18  055 


83 


Destroy  this  report  when  no  longer  needed.  Do  not 
return  it  to  the  originator. 


The  findings  in  this  report  are  not  to  be  construed  as  an 
official  Department  of  the  Army  position  unless  so 
designated  by  other  authorized  documents. 


The  contents  of  this  report  are  not  to  be  used  for 
advertising,  publication,  or  promotional  purposes. 
Citation  of  trade  names  does  not  constitute  an 
official  endorsement  or  approval  of  the  use  of  such 
commercial  products. 


security  classification  of  this  page  (Whwt  Dmtm  Entmred) 


REPORT  DOCUMENTATION  PAGE 


t.  REPORT  number 


Miscellaneous  Paper  CERC-83-2 


4.  TITLE  (mtd  Subtittm) 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


3.  RECIPIENT'S  CATALOG  NUMBER 


DEVELOPMENT  OF  A  NUMERICAL  SOLUTION  TO  THE 
TRANSPORT  EQUATION;  Report  2:  COMPUTATIONAL 
PROCEDURES 


7.  AUTHORC*7 

Richard  A.  Schmalz,  Jr. 


5.  TYPE  OF  REPORT  «  PERIOD  COVERED 


Report  2  of  a  series 


S.  PERFORMING  ORC.  REPORT  NUMBER 


8.  CONTRACT  OR  GRANT  NUMBERCaJ 


12.  report  date 

September  1983 


13.  NUMBER  OF  PAGES 


9.  PERFORMING  organization  NAME  AND  ADDRESS 

U.  S.  Army  Engineer  Waterways  Experiment  Station 
Coastal  Engineering  Research  Center 
P.O.  Box  631,  Vicksburg,  Miss.  39180 


11.  CONTROLLING  OFFICE  NAME  AND  ADDRESS 

U.  S.  Army  Engineer  District,  Mobile 
P.O.  Box  2288 

Mobile,  Ala.  36628  _ 


4.  MONITORING  AGENCY  NAME  &  ADDRESSCIf  itlHarmit  /ram  ContnlUnt  Otlict)  IS.  SECURITY  CLASS,  /of  thim  report; 

Unclassified 

ISo.  DECLASSIFICATION/ DOWNGRADING 
SCHEDULE 


1«.  DISTRIBUTION  STATEMENT  /of  Oil*  Roport; 


Approved  for  public  release;  distribution  unlimited. 


17.  DISTRIBUTION  STATEMENT  (ot  thm  9b»trmct  mttbr^d  in  Block  20,  li  dtHoront  from  Rmport) 


IB.  SUPPLEMENTARY  NOTES 

Available  from  National  Technical  Information  Service,  5285  Port  Royal  Road, 
Springfield,  Va .  22161. 


19.  KEY  WORDS  (Cantinum  on  eormroo  oido  If  nocooomry  and  Idontify  by  block  nimtbor) 

Constituent  transport  Numerical  stability,  consistency. 

Dispersion  and  convergence 

Finite  difference  approximations  Salinity 

Numerical  model  Truncation  error 


2^  ABSTRACT  roTOrwo  aMi  ff  nmitmcmofy  aoW  Idontify  by  block  niatbor) 

-This  report.  Report  2  in  a  three-report  series,  develops  the  computa¬ 
tional  procedures  for  implementing  the  numerical  approximations  to  the  transport 
equation.  The  stability  and  truncation  error  characteristics  are  presented  for 
the  schemes  selected  in  Report  1  for  the  linearized  form  of  the  transport  equa¬ 
tion.  Computational  procedures  are  next  developed  for  the  nonlinear  equation 
in  transformed  coordinates.  Numerical  approximations  near  boundaries,  the 

(Continued) 


EOfTION  OF  I  NOV  6S  IS  OBSOLETE 


_ Unclassified _ 

security  CLASSIFICATION  OF  THIS  PAGE  /Rlran  Dmim  EnlWMf) 


_ Unclassified _ 

SeCumTY  CLAMIFICATIOM  Or  THIS  PAOCf1Wl«n  DM»  _ 

20.  ABSTRACT  (Continued). 

’hydrodynamic  interface,  and  the  determination  of  dispersion  coefficients  are 
developed  for  simulating  conditions  in  Mississippi  Sound  and  adjacent  areas. 


Unclassified 


SeCUniTY  CLASSIFICATION  OF  THIS  PAGE(TWi®n  Dmim  Enfrtd) 


PREFACE 


The  computational  procedures  of  the  numerical  approximation  to  the 
transport  equation  are  reported  herein.  These  procedures  will  be  in¬ 
corporated  into  a  numerical  model  to  be  used  for  evaluating  effects  of 
proposed  dredged  material  disposal  practices  in  Mississippi  Sound  and 
adjacent  areas. 

Project  administration  and  funding  were  provided  by  the  U.  S.  Army 
Engineer  District,  Mobile.  Funding  for  publication  was  provided  by  the 
Coastal  Engineering  Research  Center  (CERC),  Waterways  Experiment  Station. 

Computational  procedures  were  developed  by  the  U.  S.  Army  Engineer 
Waterways  Experiment  Station  (WES)  in  the  Hydraulics  Laboratory  (HL) 
under  the  general  supervision  of  Messrs.  H.  B.  Simmons  and  F.  A. 
Herrmann,  Jr.,  Chief  and  Assistant  Chief,  respectively,  HL,  and  in  CERC 
under  the  general  supervision  of  Drs.  R.  W.  Whalin  and  L.  E.  Link,  Chief 
and  Assistant  Chief,  respectively  CERC,  and  Mr.  C.  E.  Chatham,  Chief, 
Wave  Dynamics  Division  (WDD) .  This  report  was  prepared  by  Dr.  R.  A. 
Schmalz,  Jr.,  WDD. 

Commanders  and  Directors  of  WES  during  this  study  and  the  prepara¬ 
tion  and  publication  of  this  report  were  COL  Nelson  P.  Conover,  CE,  and 
COL  Tilford  C.  Creel,  CE.  Technical  Director  was  Mr.  F.  R.  Brown. 


TABLE  OF  CONTENTS 


PREFACE  .  i 

LIST  OF  FIGURES .  3 

LIST  OF  TABLES .  3 

PART  I:  INTRODUCTION  4 

PART  II:  NUMERICAL  APPROXIMATIONS  FOR  THE  TRANSPORT 

EQUATION  IN  CARTESIAN  COORDINATES  .  5 

1.  Linear  Forms  of  the  Transport  Equation  .  5 

2.  Leendertse  Multioperational  Schemes: 

One-Dimensional  Analysis  .  7 

3.  Leendertse  Multioperational  Schemes: 

Two-Dimensional  Analysis  (FTCS) .  15 

4.  Leendertse  Multioperational  Schemes: 

Two-Dimensional  Analysis  (FTUS)  .  18 

5-  Spread  Time  Derivative  Schemes  .  19 

6.  Formal  Truncation  Error,  Eigenvalue,  and  Complex 

Propagation  Factor  Analysis  .  27 

7.  Flux-Corrected  Transport  .  37 

PART  III:  NUMERICAL  APPROXIMATIONS  FOR  THE  TRANSPORT 

EQUATION  IN  TRANSFORMED  COORDINATES  .  43 

1.  Development  of  the  Transformed  Equation  .  43 

2.  Leendertse  FTCS  Multioperational  Scheme  .  49 

3.  Leendertse  FTUS  Multioperational  Scheme  .  53 

4.  The  Spread  Time  Derivative  STCS  Scheme .  58 


1 


Page 


5.  Flux-Corrected  Transport  .  60 

6.  Additional  Flux-Corrected  Transport  Limiters  .  68 

PART  IV:  ADDITIONAL  NUMERICAL  CONSIDERATIONS  .  70 

1.  Approximations  Near  Solid  Boundaries  .  70 

2.  Approximations  Near  Open  Boundaries .  73 

3.  Tridiagonal  Matrix  Solution  .  77 

4.  Hydrodynamic  Interface  .  81 

5.  Dispersion  Coefficient/Determination  .  83 

REFERENCES .  86 

APPENDIX  A:  LEENDERTSE  1 -DIMENSIONAL  CASE  .  A1 

APPENDIX  B:  2-DIMENSIONAL  CASE,  WITH  Pe  =  10 .  B1 

APPENDIX  C:  2-DIMENSIONAL  CASE,  WITH  Pe  =  108  .  Cl 

APPENDIX  D:  PROGRAM  LISTING  .  D1 


2 


1  Space  staggered  finite  difference  grid  in  trans¬ 
formed  coordinates  47 

2  Datum  convention  employed  within  the  space  staggered 

grid  system  48 

3  Open  boundary  specification  of  S  75 

LIST  OF  TABLES 

Table  No.  _ Description _  Page 

I  Leendertse  Flow  Conditions  15 

II  Time  Level  n  Taylor  Series  Expansions  (Third  Order)  28 

III  Time  Level  n+I  Taylor  Series  Expansions 

(Third  Order)  29 

IV  Time  Level  n  Taylor  Series  Expansions 

(Fourth  Order)  30 

V  Time  Level  n+1  Taylor  Series  Expansions 

(Fourth  Order)  31 

VI  Truncation  Error  Analysis  Results  32 

VII  X-Sweep  Modifications  FTUS  56 

VIII  Y-Sweep  Modifications  FTUS  57 

IX  Tridiagonal  Matrix  Solution  Procedure  82 


3 


DEVELOPMENT  OF  A  NUMERICAL  SOLUTION  TO  THE  TRANSPORT  EQUATION: 

Report  2  COMPUTATIONAL  PROCEDURES 

PART  I :  INTRODUCTION 

This  report  develops  in  detail  the  numerical  approximations  to  the 
transport  equation.  In  Part  II,  the  linearized  form  of  the  equation  is 
the  subject  of  investigation.  The  stability  and  truncation  error  char¬ 
acteristics  for  the  schemes  proposed  in  the  first  series  are  developed. 
In  Part  III,  the  nonlinear  equation  in  transformed  coordinates  is  con¬ 
sidered.  The  schemes  developed  in  the  first  part  are  extended  to  the 
nonlinear  equation.  In  Part  IV,  the  numerical  approximations  near 
boundaries,  the  hydrodynamic  interface,  and  the  determination  of  dis¬ 
persion  coefficients  in  terms  of  flow  field  properties  are  developed. 

This  report  outlines  the  development  of  the  salinity  algorithm. 


The  next  step  is  the  numerical  implementation  of  these  procedures. 


PART  II:  NUMERICAL  APPROXIMATIONS  FOR  THE  TRANSPORT  EQUATION 
IN  CARTESIAN  COORDINATES 


A  Cartesian  coordinate  system  is  employed  in  all  developments  pre¬ 
sented  in  this  part.  The  stability  and  truncation  error  of  the  proposed 
numerical  approximations  are  investigated  for  the  linearized  transport 
equation.  In  this  manner  the  most  favorable  schemes  may  be  determined 
prior  to  programming  a  numerical  experimentation.  Unfortunately,  the 
transport  equation  is  nonlinear  and  no  formal  method  of  analysis  exists 
to  determine  the  appropriateness  of  numerical  schemes.  We  follow 
standard  numerical  practice  and  assume  schemes  which  possess  favorable 
computational  attributes  for  the  linearized  transport  equation  will  also 
be  suitable  for  the  nonlinear  equation. 

We  therefore  develop  linear  forms  of  the  transport  equation  fol¬ 
lowed  by  investigation  of  several  numerical  schemes  to  this  form  of  the 
equation.  The  schemes  considered  are  the  Leendertse  [1]  multioperational 
scheme  employing  forward  time  and  centered  space  derivatives  (FTCS). 

The  use  of  upwind  space  differencing  within  the  Leendertse  multiopera¬ 
tional  scheme  is  next  investigated.  The  scheme  thereby  obtained  is 
known  as  the  forward  time  upwind  space  (FTUS)  scheme.  We  next  investi¬ 
gate  several  spread  time  derivative  (STCS)  schemes  and  select  the  most 
favorable  for  further  development. 

1 .  Linear  Forms  of  the  Transport  Equation 

Let  us  consider  the  two-dimensional  depth  integrated  transport 


equation  as  follows: 


(1.1) 


(h»)  ♦  Ij  (I.US)  t 


(hvs)  =  |-  (hX^  |1)  ,  |-  (hKy  I) 


where 


h  =  total  water  depth 

u,v  =  depth  average  velocity  components  in  the  x  and  y  direc¬ 
tions,  respectively 

t  =  time 

x,y  =  Cartesian  coordinates 
s  =  constituent  concentration 

K  ,K  -  Effective  dispersion  coefficients  in  the  x  and  y  direc 
*  ^  tions,  respectively  (note  the  *  notation  has  been 

dropped) 


Equation  1.1  represents  the  conservative  form  of  the  transport  equation. 
To  derive  the  nonconservative  form,  we  expand  the  left  hand  side  of  1.1 
to  obtain  (noting  h  =  r)  "  Zj,): 


(3n  ...  3s  .  3(hu)  ...  3s  .  a(hv)  .  Ss 

)  -k-  ®  ^ »“  gj  +  -k-  *  5? 

Since  the  bottom  is  rigid,  Bz^/dt  =  0  .  Using  the  continuity  relation 
3r|/3t  +  8(hu)/3x  +  8(hv)/9y  =  0  and  collecting  terms  we  obtain 

^u)  3(hv)\  +  u  3s  +  V  (13) 

'‘\8t  8x  3y  /  \8t  Bx  ^  By)  ^ 


Then  finally  the  left  hand  side  of  Equation  1.1  becomes 


.  /8s  .  3s  .  8s \ 

^  “  3i  ^  3^) 


(1.4) 
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We  now  rewrite  the  transport  equation  for  two  important  special  cases. 
In  Case  I  we  assume  h  is  constant  and  obtain 


9s 

9t 


+ 


(1.5) 


This  result  is  also  obtained  if  PP^I^bl  ‘  Case  II  we  assume  h  , 

K  ,  and  K  are  all  constant  and  obtain 
X  y 


9s 

9t 


+ 


S.2 

9  s 


9  s 


u|^+v|^=K  ^+K  , 

9x  9y  X  3^2  y  ^^2 


(1.6) 


We  note  that  in  Equations  1.5  and  1.6,  although  u  ,  v  ,  s  are  depth 


integrated  quantities  and  and  are  effective  dispersion  coeffi¬ 


cients,  the  form  of  the  equations  hold  for  instantaneous  velocity  or  time 
averaged  (over  the  turbulence)  velocity  as  well.  In  fact  Equation  1.6 
or  its  one-dimensional  form  is  often  used  since  for  constant  velocity  it 
becomes  a  linear  equation.  Therefore,  von  Neumann  stability  analysis  may 
be  employed  to  analyze  the  characteristics  of  numerical  approximations. 


2 .  Leendertse  Mult ioperational  Schemes:  One-Dimensional  Analysis 

The  following  one-dimensional  transport  equation  is  employed  to 
determine  the  dissipative  and  dispersive  properties  of  the  multiopera- 
tional  scheme  [ 1 ] . 

+  u  -  D  =  0  (2.1) 

9t  9x  „  2 
9x 

u  ,  D  constant 

A  multioperational  analog  to  Equation  2.1  is  written  as  follows: 


-  P"  *  i;  u[(I  .  p)pj:j  -  *  CP  .  DP-;] 

■  ~2  H'l  -  ^ 


n+2  „n+l  .  At 
2Ax 


uj^d  +  a)Pj+j  -  2ap"'^^  +  (a  -  Dpj^jj 

-  B  -  <  *  '>?-!)  = »  <--> 


where 

p“  =  p(jAx,nAt) 

ot  =  -1  ,  0  ,  1  (Note  a  =  -1  for  backward  difference  in  space 

0=0  for  centered  difference  in  space 
0=1  for  forward  difference  in  space) 

The  solutions  to  the  2.2a  and  2.2b  are  expressed  by  a  Fourier  series 


P(x,t)  =  £  p*  exp  li(a  X  +  w  t)l 

1  Ul  111  III 

in=l 


(2.3) 


00 

p“  =  p(jAx,nAt)  =  P*  exp  Ii(a^jAx  +  w^nAt)) 
m=l 


where 


w  =  frequency 


a  =  wave  number 


p*  =  complex  constant  for  each  m 


i  =  fl 

Considering  only  one  general  term  in  Equation  2.3  due  to  the  linearity 
of  Equation  2.2  and  substituting  in  Equation  2.2a  we  write: 


^^^i[a(j±l)Ax+w(n+l)AtI  _  ^n+l  _  ^n+l^tioAx 
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We  thus  are  in  a  position  to  determine 
from  Equation  2.2a: 


n+l  ,  n 

j 


Therefore  we  obtain 


n  .  uAtTz-i  \  n 

Pj  ^  ^ 


n+l^ioAx  _  2^pn+l  ^  ^ ^pn+l^-iaAx j 


DAt 

(Ax)‘ 


^  pnn^-ioAx\  ^  „  ,2.4) 


which  may  be  rewritten  as  follows 


.  n+l 
AjP. 


(2.5a) 


(2.5b) 


In  order  to  simplify  \ 


recall  sin  0  =  (e^®  -  e  ^®)/2i 


therefore 


(2.6) 


Regrouping  ,  with  the  above  relations  in  mind 


A 


1 


uAt  ,  iaAx 

2A^ 


-iaAx 

e 


.  uAt  ,  iaAx  ^  -iaAx 


2) 


DAt  f  iaAx  ^  -iaAx.  ^  , 
— X  (e  -  2  +  e  )  +  1 

Ax 


(2.7) 


If  we  substitute  a  general  term  in  Equation  2.2  in  Equation  2.2b  we 
obtain 


9 


n+1  ^  Atu 

''j  ^  2S 


(1  -t  Ol)p” 


n+1  ioAx  -  n+1  .  ,  n+1  -laAx 

^  e  -  20(pj  +  (a  -  DPj  e 


DAt  /  n+1  iaAx  _  n+1  .  ^n+1  -iaAx 
— ^  (p,  e  -  2p,  +  p,  e 


)  =  0 


Which  may  be  written  as  follows: 


Vj 


Where  ,  ,  luAt  t  k  \  ^  2uAta  .  2 

A._  =  1  -  —  sin  (aAxj  +  -  sin 

2  Ax  Ax 


/a£^\  4DAt  .  2  /aAx 

{-T 


Defining  the  entire  transfer  process  as 


^n+2  _  .  ^n+1  _  2  n 
p  .  -  A,_p  .  =  ^  p  . 


We  obtain  the  amplification  factor  ^2/^^  ~  ^  which  for  stability 

(Ml  1  • 

Thus 

,  .  /zuAta  4DAt\  .  2  /aAx\  iuAt 
1  +  [-T7. - ^  sin  Hr-1  -  sin  (aAx) 


1  -  ^  sin‘  Sin  (aAx) 


We  observe  in  Equation  2.11  for  centered  space  differences 
(a  =  0)  and  |x|  <  1  .  For  backward  space  differences  (a  =  -1)  and 
u  >  0  1^1  1  1  »  while  for  u  <  0  the  cell  Peclet  number  must  obey  the 


following  relation  for 


<  1  . 


Pe  =  <  2 

c  D  - 


(2.12) 


For  forward  space  differences  (a  =  1)  and  u  <  0  1  ^1  1  1  »  while  for 


0 


Leendertse  notes  that  the  general  solution  to  Equation  2.1  may  be 


expressed  as 


p(x,t)  =  p*  exp  (i(ax  +  wt)] 


Substituting  Equation  2.14  into  Equation  2.1  we  obtain 


2  2 

iwp(x,t)  +  iuap(x,t)  -  Di  a  p(x,t)  =  0 
2 

w  +  ua  -  iDa  =  0 
w  =  o(iDa  -  u) 

and 

2 

p<'x,t)  =  p*  exp  [ia(x  -  ut)]  exp  (-Da  t) 


(2.14) 


(2.15a) 


(2.15b) 


We  observe  then  that  there  is  a  relationship  between  the  temporal 
frequency  and  the  spatial  frequency.  As  a  result,  the  complete  solution 
may  be  written  in  terms  of  the  spatial  frequency  solely.  For  a  time  pe¬ 
riod  At  ,  each  Fourier  component  is  decreased  in  amplitude  by 
2 

exp  (-Da  At)  and  is  propagated  a  distance  uAt  . 


In  the  computational  procedure  a  different  relationship  exists. 

The  eigenvalve  or  amplification  factor  may  be  used  to  study  the  dissipa¬ 
tive  and  dispersive  effect  of  the  computational  procedure  by  the  use  of 
the  concept  of  the  complex  propagation  factor. 

The  propagation  factor  is  expressed  in  terms  of  the  dimensionless 

2 

parameters  m  =  (L/Ax)  ,  D'  =  (DAt/Ax  )  ,  and  U  =  (uAt/Ax)  .  It  is  de¬ 
fined  as  the  complex  ratio  of  the  computed  wave  to  the  prototype  wave 
after  an  interval  in  which  the  prototype  wave  propagates  over  its  wave¬ 
length.  The  modulus  of  the  propagation  is  a  measure  of  the  decay  of 
the  amplitude  during  computation,  while  the  argument  is  a  measure  of  the 
computed  phase  shift. 

To  determine  the  factor  we  use  the  following  previous  results  for 
the  computed  solution.  We  consider  the  case  a  =  0  corresponding  to  the 
use  of  centered  space  differences.  We  note  from  equation  2.13  for  a  =  0 


^_l“b-ic  ,  ADAt  .  2  aAx  .  uAt  \ 

A  -  ^  ^  where  b  =  — ^  sin  — and  c  =  sin  (aAx) 

Ax'^ 


-n+2  _  ,„n 
p .  =  Ap . 

J  J 


b  =  AD*  sin^  c  =  U  sin  (aAx)  (2.17) 


From  the  solution  of  the  PDE  in  Equation  2.15  itself 


p(jAx,t  +  2At) 


=  p*  exp  |ia(jAx  -  u(t  +  2At))J  exp  ^-Da^(t  +  2At)J 
p(jAx,t  +  2At) 

=  p*  exp  ^ia(jAx  -  ut)J  exp  (-Da^t)  exp  (iau2At) 


exp  (-Da  2At)  (2.18) 


p(jAx,t  +  2At)  =  p(jAx,t)  exp  (-Da  2At)  exp  (iau2At) 
p(jAx,t  +  2At)  =  p(jAx,t)A^ 


12 


The  complex  propagation  factor  is  then  given  by  the  following  relation 


T 

m 


where 


n 


L 

m 

uAt 


Let  us  expand  Equation  2.19  using  Equations  2.17  and  2.18 


(2.1 


au2At  =  ^  (u2At) 


471  uAt  _  471 
m  Ax  ~  m 


E 


Note 


=[o ' 


C.  =  I  (1  -  4D'  sin^  +  u2 


2  25j  tan'‘  (2n/m)/ (I-AD'  sln^  (n/»))] 


C2  =  |(l  +  4D’  sin^ 


sin^  —  j  ^  (2n/m)/ (1+4D'  sin^  (n/m)^ 


1_  ('  - 


4D'  sin‘  +  U*  sin 


.  ,.„2  (n\\  -  ^  „2  .,.2  2n  i 
m 


[tan"'  U  sin  (2n/m) 

[.  1-4D'  sin^  (n/ro) 


-  tan 


*1  U  sin  (27l/m 
1+4D'  sin^  (n/ 


Ll 

/m)J 


„2  .  2  2n 

sin  — 
m 


^  (1  +  4D'  sin^ 

We  therefore  may  rewrite  Equation  2.20  in  final  form  by  defining  tempo¬ 
rary  variables  a  =  1  -  4D’  sin^  n/m  ,  b  =  1  +  4D'  sin^  n/m  ,  and 
c  =  U  sin  (27t/m)  .  Note  n  =  L^/uAt  =  mAx/uAt  =  m/U 

m/2U 


.2  a.  2 

a  +  c 

k2  a.  2 

T  =  D  c  1 

”  ^-2o/(2n/m)2 


j-  (c/a)-tan"^  (c/b)-(47t/m)uj'"^^^ 


(2.21a) 


i6 

'^'m  “  "  (2.21b) 

The  plot  of  versus  m  is  known  as  the  modulus  of  the  propagation 

factor.  The  plot  of  0^  versus  m  is  known  as  the  argument  of  the 
propagation  factor.  An  alternate  means  of  considering  T  is  given  by 
Leendertse  as  T(aL)  where  T[(2n/mAx)L]  or  T[  (27T/mAx)Ax]  =  T(27t/m) 

=  Tm  .  L  =  Ax  is  a  characteristic  length  equal  to  the  grid  size.  Al¬ 
though  we  have  not  shown  the  above  plots  here,  Leendertse  comments  that 
amplitude  and  phase  characteristics  of  the  multioperational  scheme  are 

good  for  m  >  10  ,  L  >  lOAx  . 

~  m  — 

In  the  simulation  of  Jamaica  Bay  Ax  =  500  ft*  and  L  >  5000  ft 

m  - 

*  To  convert  from  feet  to  meters,  multiply  by  0.3048. 
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For  wavelengths  less  than  5000  ft  the  amplitudes  will  be  amplified.  The 
flow  conditions  considered  for  initial  testing  are  given  in  Table  I. 


Table  I.  Leendertse  Flow  Conditions 


u  0.1  0.2  0.5 

D  0.01  ’  0.01  ’  0.01 


1.0  V  0.1  0.2 

0.01  D  0.04  0.04 


0.5 

0.04 


1.0 

0.04 


4,  4, 


4  4. 


4,  4, 


P  10  20  50  100 

e 


P  2.5  5 

e 


12.5  25 


3 .  Leendertse  Multioperational  Schemes:  Two-Dimensional  Analysis  (FTCS) 


The  following  two-dimensional  transport  equation  is  considered. 


an 

at 


+ 


V 


9Q  = 

ay 


(3.1) 


where 

n  =  constituent  of  concern 

u,v  =  constant  velocity  components  in  the  x  and  y  directions, 
respectively 

a  E  constant  dispersion  coefficient 
x,y,t  are  as  previously  defined 

Leendertse  [1]  employs  the  scheme  originally  proposed  by  Peaceman  and 
Rachford  [2]  for  diffusion  problems.  Namely,  for  the  X-sweep 

.  uAt  5^  .  «|t  .  v|t  ^  2|t 

For  the  Y-sweep 

I .  vAt  .  oAt  •  2\  n+1  1.  uAt  .  .  aAt  -  2\_n+l/2 


'y  =  <Vl,m  -  Vl,m>  ‘y  =  <Vl,m  *  Vl,. 


X  =  mAx 


y  =  My 


t  =  nAt 


If  we  eliminate  the  intermediate  level  ,  we  obtain 

/.  vAt  ,  aAt  .  2\  n  (.  .  vAt  -  oAt  2\„r 

-  —  6y  +  -^  6y  jn  _\1  +  —  6y  “  —  6y  jo 


(l  t  aft  6x  -  6x2) 


Expanding  Equation  3.3  we  obtain 


6x  +  ^  6x^^ 


(3.3) 


6y  -  ^  6y2)  (] 


'  ^  uAt  J5 
1  +  -y-  OX 


oAt  -  2\  n+1 
-y-  6x  jx\ 


=  6y2)  (i  -  fix  +  5ft  6x2)n" 


(3.4) 


Let  us  substitute  Ho  =  ^anAt^iyEAy^iPmAx  Equation  3.4  and  de- 

£,m 

fine  the  following  auxiliary  variables. 


®2  ~  4Ay 


,  _  oAt 

2  ■  772 

2Ay 


(3.5) 


If  we  employ,  the  results  of  Equation  2.6,  we  obtain  the  following 
expression  for  the  eigenvalue  \  of  the  numerical  approximation. 


6 


2 


II  +  2iaj  sin  sin 


2 

1  +  2ia2  sin  '/Ay  +  Ab^  sin  - 


We  note  then  jAj  <  1  ,  thus  the  scheme  is  unconditionally  stable. 

If  we  expand  Equation  3. A,  we  obtain  the  equivalent  two-dimensional 
difference  scheme. 


We  observe  further  that  Equation  3.7  may  be  rewritten  as  follows 


6^r)  +  vAt  6y  j 


uAt  5x 


(n"^'  ^  n") 


/  n+l 


-  aAt  6y 


-  aAt  6  X 


n+l  .  _n 


2  2 
.  uvAt  t  t  /_n+l  _n.>  auAt  ^  2^  /_n+l  ^n^ 

+  — ^ —  6x6y  (n  -  n  ) - 4 —  6y  6x  (n  -  H  ) 


6x^6y  -  n")  +  = 


where  6.  =  -  n"  .  We  note  the  underlined  terms  are  additional 


terms  required  by  the  factorization  necessary  to  obtain  the  multiopera- 
tional  scheme. 

4 .  Leendertse  Multioperational  Schemes;  Two-Dimensional  Analysis  (FTUS) 


A  forward  time  upwind  space  scheme  may  be  developed  by  considering 
the  following  general  space  derivative  in  operator  notation. 


=  6^  +  g  g  e  (-1,0,1) 

Ty  =  6y  +  g  ^ 


(4.1) 


Where 

T  ,T  s  general  first  derivative  operator 
X  y 

6x,6y  =  centered  first  derivative  operators  as  previously  defined 

2  2 

6x  ,6y  =  second  derivative  operators  as  previously  defined 

For  g  =  -1  ,  backward  space  differences  are  employed.  For  g  =  0  , 
the  previous  scheme  with  centered  space  derivatives  is  obtained.  For 
g  =  +1  ,  forward  space  differences  are  developed. 

If  we  replace  6  and  6  by  T  and  T  ,  respectively,  in 
X  y  X  y 

Equations  3. 2-3. 4  and  in  Equations  3.7  and  3.8  a  very  general  scheme  is 
obtained  equivalent  to  Leendertse 's  [1]  one-dimensional  analysis.  Corre¬ 
spondingly,  in  Equation  3.6  it  is  necessary  to  make  the  following  assign¬ 
ments  to  obtain  the  relation  in  Equation  4.3  below  for  the  eigenvalue  of 
the  general  scheme. 


2ia^  sin  PAx 


2iaj  sin  PAx 


4ajg  sin 


2 


2ia-  sin  y^y 


2ia-sin  y^y  “  4a.,g  sin 


2  yAy 

9 


(4.2) 


1  +  (^ajg  -  4bj)  sin^  -  2iaj  sin  pAx  1  ♦  (Aa^g  -  Ab^)  sil^  -*1^  -  2ia2  sin  yAy 

1  +  (Abj  -  Aajg)  sin^  ♦  2iaj  sin  pAx  1  ♦  (Ab^  -  Aa^g)  sin^  +  2ia2  sin  yAy  (^-3) 


in  which 


4a  fi  -  4b  = 

\  Ax  .2 


4a2g  -  4b2  = 


^  Ax 

CM 

/  \ 

(vAtg 

2o(At\ 

\Ay 

Ay^  / 

(4.4) 


We  observe  that  if  we  set  At  ■*  At/2  ,  a  ->•  g  ,  D  ^  a  in  Equation  2.11, 
for  a^  =  b^  =  o  ,  we  obtain  the  result  given  by  Equation  4.3.  Analogous 
to  the  one  dimension  case,  for  upwind  differencing  an  unconditionally 
stable  scheme  is  obtained  which  we  denote  as  FTUS. 

5.  Spread  Time  Derivative  Schemes 


Let  us  first  define  the  following  average  space  operators 


^£,m+l  ^  *^£,m-l 


(5.1a) 


^  (5.  lb) 

If  one  studies  the  relationship  between  Equations  3.3  and  3.6  and 

replaces  1  by  (2  +  (J  )/3  or  by  (2  +  p  )/3  ,  appropriately,  several 

X  y 

schemes  suggest  themselves.  In  each  case,  the  appropriate  time  deriva¬ 
tive  is  averaged  spatially  and  a  "spread"  in  space  time  derivative 
scheme  is  obtained.  Several  such  schemes  are  investigated  in  turn  below 
Intermediate  level  differencing 


If  we  replace  1  by  (2  +  ^t  the  intermediate  level  we 


obtain  the  following  relation. 
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If  we  combine  like  terms  in  Equation  5.5  we  obtain  the  following 
relation  in  which  the  additional  factorization  terms  (see  underlined 
terms  in  Equation  3.8)  have  been  omitted. 


(2  +  |J  \  /  n+1  n\ 

+  uAt6x  .  2 


2 


.  I  vAt6y(^L^)  -  I  aAt6y 


(5.6) 


+  3  vAt6yM  l-J - -  3  aAt6y  p 


»  0 


Opposite  inter- 

mediate  level  differencing 

If  we  replace  1  by  (2  +  lJy)/3  at  time  levels  n  and  n+1  and 
employ  standard  time  differencing  at  the  following  relation  is 

obtained 


/2  .  .  vAt  r  oAt  r  2\  /, 


n  ..  iiAt  .  aAt  r  2\„n+l 

fl  +  Ox - ^  6x  |n 


(.  uAt  .  ^  aAt  t  2\  /2  ^  vAt  r  .  aAt  .  2\^n 

=  )  (3  *  3  ■  T  ’y  *  -r  )i 


(5.7) 


TJ7  1-  .  .  .  n  anAt  iY£Ay  iSmAx  .  .  ^  .•  .  t  .u 

If  one  substitutes  n.  =  e  e  *  e  ^  into  Equation  5.7,  the 

jL  ^ni 

following  eigenvalue  is  obtained. 

^1  -  2ia^  sin  PAx  -  sin^  '^) 


A  = 


0  +  cos^Y^y  +  213^  sin  yAy  +  <ib^  sin^  (’  ^  2iaj  sin  pix  +  4bj  sin^ 


(5.8) 
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Since  A.  <  1  ,  this  scheme  is  also  unconditionally  stable.  The  scheme 
is  given  by  the  following  relationship. 


(l  *  afi  -  2ft  =  (f  »  ^  -  tfl  6y  ♦  2fl  6y2)n" 

0  ^  .  a|t  ,  oat  ^^2yn/2 


(5.9) 


Expanding  Equation  5.7  the  following  equivalent  two  dimensional  scheme 


is  obtained. 


2  ^y  vAt  t  oAt  c  2  ^  uAt  ,  ,  uAt  ^  .  uvAt^  ,  t 

-  +  ^  +  —  6y-^6y  +  —  6x  +  -^  6xPy  +  — ^  6y6x 

k 

2  2 
auAt  c  2^  aAt  <•  2  oAt  ,  2  avAt  2*  . 
- —  6y  6x  -  ^  6x  -  -g-  6x  My - ^  6x  6y 

.  a^At^  c  2,  2\_n+l 


+  6x-6y‘-\r) 


(2  vAt  j5  .  0(At  .  2  uAt  *  uAt  *  .. 

3  +  3^  -  —  6y  +  —  6y  -  _  6x  -  6xm 


(5.10) 


2  2  2  2  A 

auAt  c  2,  .  aAt  c  2  .  aAt  j.  2  avAt  r  2*  ,  a  At  x..2j.„4,n 

— ^ —  6y  6x  +  -3-  6x  +  ox  My  "  — 2^ —  oy  +  — ^ —  Ox  Oy^ 


Equation  5.10  may  be  written  neglecting  the  factorization  terms  as 
follows . 


/2  +  M 


3 — +  vAt6y 


aAt6y‘ 


/n+1  .  n\  1  9 

+  3  uAt6xMy(^^ - 2 — '  3  «^t6x  My 


n+1  .  _n 


j.o 


(5.11) 
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Advanced  level  differencing 

If  one  employs  spread  time  derivatives  at  the  most  advanced  levels 
in  each  sweep  and  and  utilizes  previous  procedures,  the 
following  eigenvalue  is  obtained  for  this  scheme. 


-  2ia^  sin  pAx  -  4bj  sin^  ~  ^**2 


+  £51^4?  +  2ia|  sin  PAx  *  4b|  sin^  ♦  £55_-1(Ay  ^  2ia^  sin  ,Ay  ♦  sin^ 


(5.12) 


Consider  a  stability  investigation  in  the  following  manner.  First,  note 
by  trigonometric  identity 

.  2  6Ax  1  -  cos  6Ax  .  .  2  vAy  1  -  cos  vAy 

sin  =  - 2 - ^ '2  ~  - 2 - 

Consider  yAy  =  (0)  =  pAx  ,  then  since  sin  (0)  =  0  and 
cos  (0)  =  1  ,  Equation  5.12  becomes 

A  =  I  I  Al  >  1 


and  the  method  is  unstable.  Therefore,  this  scheme  will  not  be  further 
considered . 

Retarded  level  differencing 

If  one  employs  spread  time  differencing  at  the  most  retarded  time 
level  in  both  sweeps,  the  eigenvalue  is  given  by  the  following  relation. 


k  = 


(?  ,  c^ax  .  Sin  PAX  -  4b,  gf^Xl  "  '^‘*2  '  ^'>2  f') 

(l  ♦  2ia  sin  yAy  +  4b  sin^  f’  +  2ia  sin  PAx  +  4b  sin^ 


(5.13) 


This  scheme  is  unconditionally  stable  and  is  given  by  the  following 
relation. 


.  uAt  t  aAt  r  2\  ,  vAt  .  o(At  .  2\  n+1 

^1  +  —  6x  -  —  6x  J  ^1  +  —  6y  -  —  6y  In 

_  ( 2  .  ^  vAt  -  .  oAt  «  2\  /2  .  *^x  uAt  .  ^  oAt  .  2\  „n 

'  \3  ♦  3  -  ^  ^  j(5  ^  3 - T  ^ 


(5.14) 


The  sweep  equations  then  become 


c  aAt  f.  2\  n+1/2  /2  ^y  vAt  t 

6x  -  ^  6x  jn  V  =  (3  ^  —  6y 


4.  .  21 

+  -Sy 


(5.15) 


L  .  vAt  *  aAt  ,  2\  n+1  _  (2  ^x  uAt  j.  .  aAt  ^  2\„n+l/2 

(1  +  —  6y  -  ^  6y  jn  =3-^3 - ^  6x  +  ^  6x  jn 


If  Equation  5.14  is  expanded  the  following  relation  is  obtained. 


I  2  2 

1  4.  ,  aAt  jt  2  ,  vAt  «  .  uvAt  *  -  avAt  ►  2* 

1  +  -y-  Ox  -  -y  Ox  +  -y  oy  +  — ^ —  OxOy  -  — ^ —  Ox  Oy 

k 

aAt  2  auAt^  *2,  .  a^At^  t  2^  2i_n+l 

-  -y  oy - 4 —  6y  5x  +  — —  6x  6y  In 


/ 4  .  2  vAt  -  ^  aAt  .2,2  ^  ^y^x  vAt  , 

=  9  5  ^  ~  ~  *  9  ^  *  9 - r 


aAt  .  2  uAt 


auAt  .  2. 


+  -y  6y  -  —  6x  -  -y  6x|jy  +  ^  6y6x - —  6y  6x 

.  aAt  t  2  ,  aAt  ,  2  ovAt^  ^  2^  .  a^At^  ,  2,  2\_n 

+  —  6x  +  -y  6x  M - ^  6x  6y  +  — y-  6x  6y  Iq 


(5.16) 
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Equation  5.16  may  be  recast  into  the  following  form  (ignoring  At 
factorization  terms). 

4  +  2(|J  +  |j  )  +  p  (j  ^n+1  (2  +  u  ) 

„n+l  ^x'^y  n  ^  n  ^  '  ^x  n 

n  -  - ^  n  +  vAt6y  +  — g - n 

n+1  (2  +  p  )  „  rt”"*"!  (2  +  M  ) 

+  uAt6x  +  - g — ^  n"  “  oAtAy  +  - g — — 


-  aAtSx  -^2 —  ■*■ 


n+1  (2  +  M„) 


n"  0 


Comnlete  time  level  differencin 


If  one  employs  spread  time  differencing  at  all  time  levels  in  both 
sweeps,  the  scheme  eigenvalue  is  given  by  the  following  relation 

(I  *  .  2ia  sin  PAx  -  4b  sin^  (I  ♦  Y^y  .  2ia,  sin  vAy  -  4b,  sin^  ^ 


t  y  +  2ia2  sin  y^y  +  4b2  (3  3  ^  +  2iaj  sin  pAx  +  4bj  sin 

This  scheme  is  unconditionally  stable.  Corresponding  to  Equation  5.18, 
the  scheme  becomes 


(2  .  ^x  uAt  -  ^  oAt  s,  2\  (1  ^  vAt  .  .  0(At  -  2\„n 

(3*3 - r  TT j(3  *  3  ■  T  ■*  T 

/2  ^  ^  vAt  ,  aAt  .  2\  /2  ^  uAt  .  aA£  .  2\  nH 

=  (3  *  3  *  —  jis  *  r  *  T 

In  multioperational  form,  the  scheme  is  given  by 


(5.19; 


(I  4  ,  -  ay  -  Sft  =  (i  4 


uAt  jj  aAt  ►  2) 
6x  +  -^  6x 


/2  ^x  uAt  -  aAt  2\  n+1  /2  .  vAt  ^  .  aAt  ^  2\_n+l/2 


two  schemes  employ  spread  time  derivatives  in  both  directions.  For  a 
two-dimensional  computation,  the  first  two  schemes  appear  less  desirable 
than  the  last  two  schemes.  These  last  two  schemes  are  therefore  further 
investigated  within  a  formal  truncation  error  analysis. 

6 .  Formal  Truncation  Error,  Eigenvalue,  and  Complex  Propagation 
Factor  Analysis 

In  order  to  compare  the  schemes  developed  with  respect  to  trunca¬ 
tion  error,  Taylor  series  expansions  were  developed  for  the  constituent 
terms  common  to  all  schemes.  In  Tables  II  and  III  the  expansions  are 
carried  through  third  order,  while  in  Tables  IV  and  V  the  expansions  are 
carried  through  fourth  order  terms.  Substituting  the  appropriate  expan¬ 
sions  for  the  terms  in  each  scheme,  it  is  shown  that  all  schemes  are 
consistent  with  the  linearized  transport  equation.  The  order  of  the 
principal  truncation  error  is  given  in  Table  VI  for  each  scheme. 

We  note  that  the  complete  time  level  differencing  spread  time 
derivative  scheme  is  truly  second  order.  Therefore,  it  is  the  more  ac¬ 
curate  of  the  two  spread  time  derivative  schemes  and  will  be  the  subject 
of  further  numerical  development. 

The  Leendertse  multioperational  schemes  in  tandem  form  a  lower 
order  (FTUS)  and  higher  order  (FTCS)  pair,  which  may  be  developed  within 


flux  corrected  transport. 


9  A 


riil:- 


T 


Table  III.  Time  Level  n+1  Taylor  Series  Expansions 

(Third  Order) 


n+1 

n 

'  £,m 

n+1 

n 

= 

■'£,1 

r  n+1 

1 

^  Ho  = 

y  '£,m 

2Ay 

n+1 

n 

»^x^£.m  = 

^'£. 

6  = 

1 

2  2  2  2 
f  At  3  r 

'  2  2'  2 
3y  ^  at 


2  3 

^  At  a-^n 

2  2 

ay'^at 


9n  .  Ay^  a^n  .  .  ..2  a^n 


2  _2  .^2  _2 

i_  a  n  .  At  a  r 

'  2  2'  i 

ax  ^  at 


2  2  3'  3 

ax"^at  at 


X  £,m  2Ax  \  ax 


2Ax  ^  t  2AxAt  ^  ^  +  AxAt^ 

2Zix  i-  2Axat  +  33  +  ^^t 

k  ax  axat  , 


M  u  n”*'  =  n"  t  At  5D  +  sffl  ^  afn  ^  ^  3^ 

'•xVj,,  Atg^  2,  21  3^2  2!  ^^2  3! 


2  .3. 


.2  ^3. 


1  i 

^y^x'''2,m  2Ay\ 


i?  ’ ' 


=  -L/^ 

,  m  2  Ax  r 


2Ax  +  2AxAt  ^  2  +  AxAt 

ax  axat 


2  2 

atax  ^ 

atay' 

2  a^n  .  .  2. 

- 1-  +  Ax  Ay 

’  0 

ayat 

Bx  dy 

2  a^n  .  .  *  2  a^n 

- h;  +  AxAy  - ^ 

axat 

dxdy 

3  „  3 

ay  ! 


6x6yr| 


n+1  _ _ 1_ 

£,m  4AxAy 
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Ta.jle  V.  Time  Level  n+1  Taylor  Series  Expansions  (Fourth  Order) 


Table  VI.  Truncation  Error  Analysis  Results 


Order  of 

_ Scheme _  Equation  Truncation  Error* 

Multioperational 

Leendertse 


FTCS 

3.7 

0 

(At  , 

FTUS 

3 • 7  and  4 . 1 

0 

(At^, 

Spread  Time  Derivative 

Complete  Level 

5.21 

0 

(At^, 

Retarded  Level 

5.16 

0 

u. 

2  2 
Ax  ,  Ay  ) 

Ax,  Ay) 


2  2 
Ax  ,  Ay  ) 

A  2  ,  2\ 

Ax  Ay  \ 
At’  At/ 


L  A  A  *^2  .  lim  (l  -  L  I  ^  V  u  a 

\0  ASj  ,  Aa2  !♦-»  '  o  -  /  H.Aa. 

'  '  Aa .  o  ^  ^ 


where  the 


Aa . 
1 


i  =  l 


H. 

1 


are  bounded  and  L 


i  =  1 , . . .  k 

is  the  finite  difference  operator  and  L 


the  differential  operator.  2  2 

**  Ax.  Ay.  At  ^  o  -  ^  ^  = 

^  Ax,  At  -»  0  At  Ay,  At  >  0  At 

dition  does  not  hold  the  scheme  is  not  valid. 


o  . 


If  this  con 


Eigenvalue  analysis 

To  facilitate  this  analysis,  the  following  dimensionless  quantities 
are  defined.  In  the  general  two-dimensional  case,  the  dispersion  coeffi¬ 
cients  are  different  in  each  coordinate  direction  and  are  represented  by 
and  ,  respectively. 

In  previous  paragraphs,  D  and  a  have  been  used  to  represent, 
a  constant  dispersion  coefficient  in  the  one-  and  two-dimensional  cases. 


It.  _ 
nAx 

2tt 

L 

n 

2tt 

m  mAy 

U  = 

Ax 

V  = 

Ay 

9 

4n“ 

- 1- 

-f 

_2  4tt^ 

D  At 

ni  -  _ 

D  At 

n'  -  y 

2.2 
n  ,AX 

v- 

J 

%  2.2 
m  Ay 

u  -  r 

^  A  2 

Ax 

~  2 
^  Ay 

General  two-dimensional  scheme 

Consider  equation  (4.3)  which  determines  the  eigenvalue  for  the 
following  three  schemes  if  u,v  >  0  : 

(1)  g  =  -1  :  Upwind  differencing  (FTUS) 

(ii)  g  =  0  :  Centered  differencing  (FTCS) 

(iii)  g  =  +1  :  Forward  differencing  (FTFS) 


Noting  =  U/4  , 


V/4  ,  =  DV2  ,  and 


following  expression  is  obtained.  (6  =  a  ,7=0) 

n  m 


1  +  (gU  -  2D^)  sin^ 

-  i  1  sin 

^2ti\ 

1  +  (2D^  -  gU)  sin^ 

^  .  U  . 

+  1  y  sin 

1  +  (sV  -  2D’)  sln^  ‘  I  si” 


'  -  •  ^(?) 


1  +  (2D  -  gV)  sin 


.  .  V  .  /2Tr> 
+  1  —  sin 


sx 


D'/2  the 

y 


(6.2) 


SV 


Spread  time  derivative  scheme 


The  complete  spread  time  derivative  scheme  eigenvalue  given  in 
(5.18)  may  be  written  as  follows  in  terms  of  dimensionless  quantities. 


2  ^  J  u  AN  .  (A 

_ 3  _ 2  \  m  / _ X _ \  m/_ 

2  ,  w)  ,  .  U  .  /2A  ......  .  2  /  tt\ 

-z  +  - -  +  1  V  sin  ( — )  +  2D  sin  I  —  ) 

3  3  2  \m/  X  \m/ 


(^)  . 


V  .  /'2ti\  ,  .  2  /  it' 

-r  sin  I — I  -  2D'  sin  I  — 

_2 _ \n/  y _ \  n  , 

1  sin  m  *  2D'  si„2  fi' 

2  \n/  y  Vn> 


X  •  X 
sx  sy 


Computation 


The  eigenvalue  of  each  scheme,  ,  has  been  expressed  in  terms 

of  the  X  and  y  wave  numbers,  n  and  m  respectively.  The  eigen¬ 
values  are  compute  for  n  and  m  values  from  2-9  over  three  log 

cycles.  Note  L  =  nAx  and  L  =  mAy  ,  such  that  the  wavelengths  are 
n  m 

expressed  in  terms  of  the  grid  spacing  interval. 


Complex  propagation  factor  analysis 


In  order  to  compute  this  quantity,  it  is  fir-',  necessary  to  deter 
mine  the  eigenvalue  of  the  analytical  solution,  A  .  Consider 


1 


in 


Then  n  C  exp  [iS  t  +  i(a  x  +  o  y)  J  is  attempted  as  a  solution 
n  n  n  m 

Calculating  the  derivatives: 


^  =  ie  n 

3t  n 


2  2 

(6  +UO  +vo)i=-oD  -oD 
n  n  m  n  x  m  y 


3  =  ifo^D  +  o^D  )  -  uo  - 

n  \  n  X  my/  n 


6  =o(iDo  -u)+o(iD0  -v) 

n  n  X  n  m  y  m 


Thus  the  solution  may  be  expressed  as  follows. 


C  "  C  exp  S  J  -  i(uo  +  vo  )  t  +  i(o  x  +  a  y)> 

n  |L\nx  my/  n  mj  n  m  l 

C  'V  C  exp  -(o^D  +  a^D  V  exp  [io  (x  -  ut)  +  io  (y  -  vt)  ] 

n  ^|_\nx  my/J  ^  n'  m 


C  v  c  exp 
n 


exp  [ia^(x  -  ut)] 


(-o^D  t^  exp  [io  (y  -  vt)] 


Therefore  the  analytical  eigenvalue  is  given 


=  exp  ^-a^D^At^  exp  (-ia^uAt)  exp  ^^-a^D^At^  exp  (-ia^vAt) 


(6.9) 


X 

a 


In  terms  of  dimensionless  quantities,  we  finally  obtain: 
|eKp  [-(^)  dJ  exp  (-1  ^  u)  I 

I  '’‘f  [-(^)  “ ■  >.x  •  *e 


(6.10) 


The  complex  propagation  factor,  C  ,  is  computed  as  follows. 


(6.11) 


where  M  =  ^  and  N  =  -^ 

We  note:  uMAt  =  L  =  mAx 
m 

vNAt  =  L  =  nAy 
n 

A  computer  program  has  been  written  to  determine  both  the  eigenvalue 
and  complex  propagation  factor  for  the  previous  schemes  at  different 
values  of  x  and  y  wave  numbers,  m  and  n  ,  respectively. 

Initially  the  one-dimensional  case  with  U=0.2,  V=0.0, 

=0.01  ,  and  =  0  was  computed  and  results  compared  with 
Leendertse's  analysis.  The  results  matched  exactly. 
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Next,  a  two  dimensional  case  with  U  =  1.0  =  V  and  D  =  D  =  0. 

X  y 

was  considered.  Finally,  the  following  prototype  condition  case  was 
considered . 


,,  _  uM  _  (3  fps)(360  sec)  _  „  _  vM  _ 

Ax  ~  (4000  ft)  "  "  Ay  ■ 


0.  =  =  (!00_ft^/sec)(360  s^e)  ^ 

^  Ax  (4000  ft)"^  Ay^  ^ 


(6.12 


Results  for  the  three  cases  above  and  the  program  listing  are 
presented  in  Appendices  A-D. 


7.  Flux-Corrected  Transport 


In  the  implementation  of  this  method,  both  higher  and  lower  order 
in  space  schemes  are  considered.  The  schemes  are  written  in  the  follow 
ing  flux  formats. 


nj  =  n"  -  (AxAy)“^(Fj 

d  ,m  ,m  \  a 


^  -  F^  +  F^ 

)l+l/2,m  Z-l/2,m  ll,m+l/2 


where  t  =  nAt  ,  x  =  mAx  ,  y  =  ilAy 

=  concentration  at  location  ()l,m)  at  time  level  n 

ic  ,m 

Ax  =  X  space  step 
Ay  =  y  space  step 


I  =  general  index  at  time  level  n+1  which  we  set  equal  to 
H  and  L  for  the  higher  and  lower  order  schemes, 
respectively . 


^4+1/2, m+;l/2 


fluxes  through  the  appropriate  cell  faces  of  cell 
(£,m).  Form  dependent  upon  the  finite  difference 
formulation. 


We  observe  from  (7.1)  that  the  difference  between  the  higher  and  lower 
order  scheme  at  (£,m)  may  be  written  as  follows: 


=  -(AxAy) 


H 

£+1/2, m 


F 

£+1/2, m 


) 


/f”  -  U  /'f” 

V  £-l/2,m  £-1/2, m/  \  £, 


m+1/2  £,in+l/2> 


(7.2) 


■(f« 

V  ^.m- 


1/2 


£  ,m-l/2 


)] 


Note  this  difference  is  expressed  as  an  array  of  fluxes  between  adjacent 

grid  points  and  is  the  condition  required  to  implement  flux-corrected 

transport.  We  next  develop  the  expressions  for  the  above  fluxes  for  the 
H  L 

higher  (F  )  and  lower  (F  )  order  schemes. 

For  the  higher  order  scheme  we  employ  the  FICS  scheme  written  be¬ 
low  in  which  the  factorization  terms  necessary  for  the  multioperat ional 
method  are  underlined. 

n+1  n  vAt  ,  /  n+1  .  n\  uAt  .  /  n+1  .  n\ 

n„  -  n  -  —  +  n  )  -  —  «x(l„  +  n  ) 

(7.3) 

2  2 
uvAt  r  ,  (  n+1  n\  ,  auAt  „  2,  /  n+1  n^ 

-  —r~  -  h  j  ^  6y  6x(^n^^  -  n  ) 

.  avAt^  ,,2,  /  n+1  n\  a^At^  „  2.  2/  n+1  n\ 

+  — ^  6x  6y  -  n  ) - —  6y  6x  -  n  j 
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Ignoring  the  factorization  terms  (7.3)  may  be  written  in  the  form  of 
(7.1).  The  total  fluxes  are  presented  as  the  sum  of  the  advective  and 
diffusive  fluxes  by  defining: 


pH  _  HA  HO 

8.+l/2,m  il+l/2,m  1+1/2, m 


pH  _  pHA  HO 

^,m+l/2  "  «,,m+l/2  £,m+l/2 


(7.4) 


where , 


^)i+l/2  m+1/2  "  higher  order  scheme  fluxes 
HA 

^?.+l/2  m+1/2  ^  higher  order  scheme  advective  fluxes 
HO 

F^+j^/2  ^  higher  order  scheme  diffusive  fluxes 

Expanding  Equation  (7.3)  using  the  definitions  listed  at  the  top  of 


page  16  one  then  obtains: 


F.  , ,  =  vAtAx  — 

<l+l/2,m  I 


4, m+1/2 


=  uAtAy 


/  n+1  ,  n\  /  n+1  ,  n\  1  / 

2  2 

7  n+1  .  n\  /  n+1  ,  n\  / 

" "  A,.  L 

2  2 


(7.5) 


(7.6) 


Sl+l/2,m 


-aAtAx(n^  ,  +  n?..  _  “  ~  n?  )/2Ay  (7.7) 


t,HO  .  .  /  H  .  n  H  n  \ 

F.  =  -aAtAxIn.  +  n.  "  -  n.  ,  j/2Ay  (7.8) 

£-1/2, m  \  £,m  £,m  £-l,m  £-l,m// 
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pHO 


-aAtAyfn^  ,  +  n?  )/2Ax 

\  H.m+l  Z,m  Z,m/I 

-aAtAy(n^  1  “  n"  ^  / It-x. 

\  il,m  S,,m  Ijin-l 


(7.9) 


(7.10) 


Next  consider  the  FTUS  lower  order  scheme  given  below.  For 
g  =  +1  ,  u,v  <  0  ,  respectively.  Factorization  terms  are  again 
underlined. 


(7.11) 


If,  as  in  the  previous  case,  the  factorization  terms  are  ignored,  (7.11) 
may  be  written  in  the  form  of  (7.1).  Total  fluxes  are,  as  before, 
presented  as  the  sum  of  advective  and  diffusive  fluxes.  Thus 


p*^ 

£+1/2, m 
F^ 

£ , m+1 / 2 


pLA  .  pLO 

£+1/2, m  £+1/2, m 


pLA  .  pLO 

£,m+l/2  £,m+l/2 


(7.12) 


ur 


where 


/T  -  lower  order  scheme  fluxes 
£+1/2  ,in+l/2 

LA 

/o  .1/0  =  lower  order  scheme  advective  fluxes 
?,+l/2,ni+l/2 

/o  .1  /o  =  lower  order  scheme  diffusive  fluxes 
i+1 / 2 ,  n^l / 2 

Expanding  Equation  (7.11)  one  then  obtains: 


7,,m+l/2 


/  n+1  ^  n 
uAtAy\ - 2'—- 


/  n-^1  ^  n\ 

uAtAy\ - X - / 

\  /£,m+l 


u  <  0 


(7.13) 


£,m-l/2 


/  n+1  ^  n\ 

(\  I 

uAtAy\ - X - }  u  >  0 

V  ^  7«,m-l 

/  n+1  ^  „n\ 

uAtAyV - x - f  u  <  0 


(7.14) 


pLA 

1+1/2 ,m 


vAt  Ax 


vAt  Ax 


n+1  ,  n 
\  +  n 


n+1  n 

\  +  n 


V  >  0 


V  <  0 


1+1,  m 


(7.15) 


£-1/2,1 


vAtAx 


vAtAx 


n+1  ,  „n 
\  +  n 


n+1  ,  n 


£-l,m 


V  >  0 


V  <  0 


(7.16) 
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PART  HI:  NUMERICAL  APPROXIMATIONS  FOR  THE  TRANSPORT 
EQUATION  IN  TRANSFORMED  COORDINATES 


The  transport  equation  is  transformed  from  x-y  space  to  -  Ct^ 
space  by  means  of  an  exponential  stretch.  Subsequently,  the  extensions 
of  the  numerical  approximations  to  the  nonlinear  transformed  transport 
equation  are  presented.  It  is  instructive  to  note,  that  even  the  lin¬ 
earized  transport  equation  becomes  nonlinear  in  transformed  coordinates. 


1 .  Deve  1  opment  of  tlie  TranFormed  Equation 


The  following  coordinate  transformation  is  considered  by  Butler 


x  =  a,  t  bjOj 


,*  ■  “A 


(1.1) 


y  =  a2  +  b20f2 


/y  - 

V  ”2  / 


(1.2) 


Then  for  an  arbitrary  hydrodynamic  variable  p(x,y,t) 


da,  ,,  r,  da„ 

3p  _  9p  _ 1  ^  _  9p  _ 2 

9x  ~  9a ^  dx  9y  9a2  dy 


9  p  ^  9 

^  2  9a, 

9x  1 


\9x/  dx 


^  /'3e_  ^ 

9a j  l9aj  dx  I  dx 


(1.3) 


^  ^  ^  ^  (a 

dx  dx  9a  j  9a  ^  ydx  ^ 


O  -  _3_  _ 

9y2  \dy}  dy 


(1.4) 

.  JL  /SiL  ^12  -  f?!e  ^2  a  ^ 

V9yy  dy  9a2  y9a2  dy  f  dy  dy  ^^2  dy  9a2  9a2  '^dy  J 


If  we  introduce  =  dx/da^  and  =  dy/da^  then 


^  9£_ 

3x  ■  Mj  3«i  9y  ^2  9«2 


^_1_  i_9^^§2 _ 9_/jL 

3x2  -  Ml  Ml  a„2  30^  3a^ 


^_i_  i_^^3e _ ^(i_ 

3y2  ^^2  *^2  3a2  ^“2  ^“2  \^2 


Considering  Equation  1,4a  in  an  alternate  manner 


2 

^  ^  _  d_  (dp  \  ^  ^  ^ 

3^2  -  3x  I  30^  dx  y  -  3x  ^30^  dx  3o,  ^,^2 


Noting  3/3x  =  (3/3aj)(daj/dx)  =  0/3aj)(l/Mi) 


2  2 
3  P  _  9  P 

2  2 
3x  3aj 


/^:iy , 

I  ao, 


Employing  previous  notation,  Equation  1.9  is  rewritten  as  follows; 


9!fi  =  ^  {L\^  +  9£_ 
3x2  3^2  (^mJ  9a^ 


:  i  fe) 


Note,  however,  from  the  relation  between  3/3x  and  3/30(j  we  obtain 


^  +  92_  ^  /1_\  L 

_  2  »  2  lu./  30^  da^  I  Mi  /  u. 
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This  relation  is  equivalent  to  Equation  !.ii. 

If  we  consider  a  hydrodynamic  variable  and  let 

be  defined  such  that 


P,*  i'.v  “  P(i*Aa„,j*Aa  ,nAt) 


(1.12) 


Theti  let  i  ,  j  ,  n  be  such  that 


n  ^2  1 

~  P  ^2  b^d^Aa^)  ,  +  b^(j*Aaj)  ,  nAt  (l.i: 

Wo  onip  I  oy  uniform  spacing  in  “  ^2  ^Pace  and  irregular  spacing  in 

x-y  space.  We  may  evaluate  the  derivatives  with  respect  to  x  and  y 
.IS  1  <1  1  1  OWS  . 


§£ 

"  =  §fi- 

n  da. 

8x 

.  .  8a, 

dx 

1 

i’'S  j" 

(1.14) 


she  re 


da,  _  I  /x  -  a, 


dx  c,b,  I  b, 


(1-C,)/Cj 


=  f(x) 


/  c\  (1-c) 

f  a,  +  b,a,  =  -4-  a, 

\*  >  V  'i'’i  ' 


^  ST  ....  = 

j" 


where 


ex  n 

8£ 

-  §£_ 

n  da2 

8y  .  . 

1 .  J 

8a2 

dy 

i-’S  j-' 

da2  _  , 

V  (1-c  )/c 

/y  -  a 

dy  C2b2 

1  '■a 

] 

(1.15) 


=  8(y) 


cA  (1-c  )  da 

g|a2  +  j  =  ^  «2  =  W  i*  " 


For  the  second  derivative  term  we  obtain 


ax^ 


1 

"  ^ 

n 

d 

1 

dx 

.  ] 

laa. 

dx 

aa. 

da 

\dx  J 

....  ( 

j  1 

1  1 
s. 

J  1 

i*,j* 

J'-) 

(1.16) 


where 


4  ^  ‘-.“D]  =  ^  ^r’  =  -c*.) 


/c 


^  - 


dOj  \< 


.dx/  j.. 


h(j*Aap 


Simila 


rly,  for 


The  underlined  terms  in  Equations  1.6  and 


ay‘ 


l.J 


1.7,  although  they  may  be  computed  exactly,  are  approximated  using  finite 
differencing  on  and  • 

Transforming  Equation  1.1  in  Part  1  in  x-y  space  to  Uj-a2 
space  we  obtain  the  following  result. 


(dus)  (dvs) 
a  a. 

(ds)  +  -  +  - - 

t  Mj  Mo 


1 

(s) 

a. 

dK  - - 

(s) 

a 

2 

dK  - - 

^^l 

- 1 

55 

_ 1 

M2 

“2  ^2 

(1.17 


where  d  is  introduced  as  the  depth  in  place  of  h 


(  =  a/at 


(  )„  =  a/aa 

1 


(  =  a/aa^ 


Equation  1.17  is  the  relation  that  is  the  subject  of  numerical  approxi 
mation.  Let  us  consider  the  space  staggered  grid  shown  in  Figure  1. 
The  datum  convention  is  illustrated  in  Figure  2. 
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Let  us  introduce  the  following  notation  as  a  prelude  to  the  approx¬ 


imations.  Define  for  an  arbitrary  variable  F  ,  where  t  =  kAt  , 

n,m 


y  =  nAy  ,  x  =  mAx  : 


-  F^ 

t  \  n,m/  n,m  n,m 

-  F*^ 

t  \  n,m/  n,m  n,m 

5  (f^  ^  =  F*^  -  F^ 

a^Vn,m/  n,m+l/2  n,m-l/2 

6  (^F*"  'l  =  F*"  -  F^ 

a^V  n,m/  n+l/2,m  n-l/2,m 

*^1  { +  F^  ^ 

— =  _  \  n,nH-l/2  n,m-l/2/ 


nH-1/2  ^n,m-l/2 


(1.18a) 


(1.18b) 


(1.18c) 


(1.18d) 


(1.18e) 


h  +  F^  '' 

—  _  V  nl-l/2.m  n-l/2.m/ 


(l.l8f) 


2 .  Leeniertse  FTCS  Multioperational  Scheme 

The  following  finite  difference  equation  is  considered  as  an  ap¬ 
proximation  to  the  nonlinear  transport  equation  (1.17) 


',’Nds)  +  6 

t  2Aa^(y^)^ 


/a  a  a  a  \ 

I  ^k+1— k+1  k+1  ^  ^k-i k  k 
\  a  s  u  +  dsu/ 


+  6 
2Aa„(u  )  a„ 
Z  Z  n  Z 


/«  a  a  a  \ 

l-4k+l— k+1  k+1  .  -fk— -k  k) 
\d  s  V  +  dsv/ 


a  6  (s’^'^b  a  6  (s*^) 

- 6  -7-T-  -tV- 

2(Aa,)^(y,)  “l  “l  ^^l^m  1 

1  1  m  L  -J 

.  ,  k+1.  ^  ,  k. 

6  (s  )  6  (s  ) 

- - 6  4V- 

2(Aa.,)^(u„)  “2  “2  ^^2^n  “2  ^^2^n 

2  2  n  *- 


(2.1) 


=  0  at  (n,m) 
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The  above  equation  is  assumed  to  be  contained  within  the  following 

multioperational  difference  equations.  For  the  linear  case  obtained  for 

(u„)  =  (y,)  =  1  ,  K  ,  K  constants  in  space  and  time,  and 

2  n  1  m  “2  ^'l 

u  ,  V  ,  d  constant  in  space  and  time,  the  constituent  intermediate 
time  level  may  be  eliminated  in  the  multioperational  approach  and  the 
total  difference  equation  obtained  equals  the  above  difference  equation 
plus  some  higher  order  in  time  factorization  terms.  The  total  differ¬ 
ence  equation  is  consistent  with  the  linear  transport  equation.  For 
the  nonlinear  case  considered,  it  is  not  possible  to  eliminate  the  con¬ 
stituent  intermediate  time  level.  Thus  the  exact  form  of  the  factoriza¬ 
tion  terms  may  not  be  determined.  However,  their  numerical  effect  may 
be  tested. 

The  approximations  for  the  X-Sweep  may  now  be  written  as  follows 


At  6 


a 


“1  “1 


6!^(ds)  +-5X  ,  ^ 

t  2Aa  (p^) 

i  t  m 


1_  l^k+1/2*  gk+l/2*^k+l/2*j 


At  6 


2Aa. (y,)^ 
1  1  m 


-Jk+1/2*  k+1/2* 

a  K 

"1 


i  m 


+  ^ 
2(y2)j^Aa2 


!2 

.  jk  k  k ; 
Vd  s  V 


(2.2) 


At  6 


2Aa2^V2^n 


d  K 


6  (s*^) 

“2 


“2  ^^2^n 


=  0  at  (n,m) 
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Expand  (2.2)  using  (1.18)  and  collect  all  terms  at  time  level  k+1/2* 


to  obtain:  (K  H  K  ) 

X 


k+l/2*  at 

2aa,(u,) 

1  In 


/  k+i/2*  .  ^  k+ 

In  ,“11  ,  +  n 

\  n,m-l _ n,tn-l  n, 


/  k+1/2*  . 

yn.nrt-l _ n. 


pH-I  n,m _ n,m/  k+1/2* 

2.  ^n,o+l/2 


/  k+1/2*  k+l/2*\ 

\ , ffi+1  ^n.m  / 


n,m/  k+i/2* 
^n, 0-1/2 


/  k+1/2*  k+1/2*' 

*  vn,m-l  ^n»m 


at  [(^n 


2aar(i.,)„ 

1  1  nt 


k+1/2* 

^n,nrt-l  n.m+1  _ j 


k+1/2*  \  /  k+1/2*  _  k+l/2*\ 

'n.m _ n.m/  V^n.nr+l _ n.a  /  j^k+1/2* 


^''l'a+1/2  ’‘n,m+l/2 


/  k+1/2*  .  .  k+1/2*  \  /k+1/2*  _  ^k+l/2*\ 

-  V'^n.m-l  ~  n.m-1  ^n.m _ n,a/  V  n,m _ n,m-l  /  jjk+1/2* 

‘“1^111-1/2  ’‘n,i»-l/2_ 


Collecting  all  terms  in  (2.2)  at  time  level  k  denoting  the  result  as 
B  ,  we  obtain:  (K  =  K  ) 

m  y 


k  _ at _ \''n+l.m  ”n+l.m  ''n.m  n.m/  k  V  n+l.m  n,ro/ 

m*  **  n.m  2a',^(p-)  2.  n+l/2,m  2. 

2  2  n  L 

/k  .  .k  \  /k  .  ^  \ 

vn-l,m  ”  ”n-l.m  ^n,(n  n.m/  k  ^n,{n/ 

2,  Vl/2.m  2. 


2(u2)n(^«2^ 


At _ I  \  n+l,n>  ^n+l,m 


+  n  -h  ifs^,  -s  I 
n,pi  n,m/  \  n+l,o  n,m/ 


^''2^n+l/2  ^n+l/2,ra 


/  k  ^  k  \  /  k 

(n,  -h,  +n  -n  )(s 

\  n-l,m  n-l,m  n,m  n,p>/  \  n 


\  n-l,m  n-l,m  n,m  n,p>/  \  n.m  n-l,m/  j^k 

*‘'2^n-l/2  ’'n-l/2,m 

In  (2.3)  we  define  -a  ,  a  ,  and  a  as  follows 

n,m-i  n,nri-l  n.m 


M 


n  ,m-l 


..\-^)k+l/2*  k+1/2*  -k+1/2* 

^  ^n,m-l/2  %,m-l/2  ^^x^n,m-l/2 

2Aai(Ui)^  2. 


AA^r+1/2*  k+1/2* 

^  ^'^n.m+1/2  %.m+l/2 


n,m+l  2Aaj^(p^)^ 


(K 

x' n,m+l/2 
^“l(^'l^m+l/2 


^  ^  ^k+1/2*  ^  At 

n,m  n,m  2Aoi^(y2^)^ 


(-|  )k+l/2*  (-i  ) 

2 


2Aa  (u  ) 

1  1  m 


k+l/2* 

^  x/n,iirfl/2 
^^l^m+1/2 


k+l/2* 

n,in-l/2 

2 


dK 

x/n,m-l/2 

^^l^m-1/2 


(2.7) 


Collecting  all  results  we  obtain  the  following  interior  equation 


for  the  X-Sweep 


a  ,  *  a  +  a  -  B  (2.8) 

n,m-l  n,m-l  n,m  n,m  n,nri-l  n,in+l  m 


The  approximations  for  the  Y-Sweep  may  now  be  written  as  follows: 


( 


X  i  /  “2  Jll.  N  At  6„  o-  «„ 

“2  V  .k+1  k+1  k+1/  2  -7k+l  „k+l  “2 _ 


(2.9) 


«  «  (  'j  At  «  f  °i 

“l  I  .k+l/2»  k+l/2*  k+l/2V.  1  >+1/2*  „k+l/2*  Jh _ 

*2AV^V  A  .  /  \  (p,). 


•  0  at  (n.a) 


Expanding  (2.9)  by  employing  (1.18)  and  collecting  terms  at  time  level 
k+1  on  the  left  hand  side  and  leaving  terms  at  time  level  k+l/2*  on 
the  right  hand  side  the  following  interior  equation  for  the  Y-Sweep  is 
obtained. 


k+1  k+1  k+1  „ 

a.  s,  +a  s  +a.T  s,,  =B  (2.10) 

n-l,m  n-l,m  n,m  n,m  n+l,m  n+l,m  n  '  ' 


where  (K  s  K  ,  K  H  K  ) 
X  ’  y 
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B  .  — (~^r7 (^pj  K+i/2*  fc+i/2* 

n  n.m  2lu^}^d,a^  |_'  /  n,»fl/2  “n,BfI/2  '  '  i»,»-l/2  “n.B-l/zJ 

I7“l  \  -  >+l/2*\  /“,  \  /  k+1/2*  .k+l/2*\1 

♦  - ^ - -  ('dlS  V  n,>H  *n,» - 1  _  j-jij  |k+l/2»  V^nj| _  n,m-l  /[ 

L'  ’■'"'"*1'*  <''l>^l/2  '  x/n,.-l/2  J 


.  Leendertse  FTUS  Multloperational  Scheme 

The  following  finite  difference  equation  is  considered  as  an 


pproximation  to  the  nonlinear  transport  equation  (1.17): 


6'^(ds)  + 


2Aa 


At  .  -^+1-  k+1  .  -4k-  k 

—7 — ^  6\d  s,  u  +  ds,  u/ 

I'nV  “A  ^  W 


At 


k+1 


ry  ».■*.■  J.  Ot~  IC 

.  _  f  1  — 7k+l  —  k+1  .  — 7k—  k , 

+  777 - 7 — ^  6  \  d  So  V  +  d  s„v 

“2\  2  2 


At 


2(Aa.)2(u,)„  “1 

1  i  m 


a 


X  /  k+1 .. 

1  ^  “1 

^+1  1  ___ 

1  1  m  1  1  m 


(3.1) 


At 


2(4<.2)'(‘‘2>„ 


“9  ^r,  “9 

-4k+l  „k+l  “2  .  -4k  ,.k  “2 

a  K  - 7 - r -  +  a  K.  — 7 - r - 

“2  ^’^2>n  “2  ^^2>n 


=  0  at  (n,m) 


The  following  upwind  difference  operators  used  in  the  above  equa¬ 
tion  are  defined  at  (n,in)  as  follows: 


“n,in-l/2 

^k 

“n,in+l/2 


f^  >0 
n.m  — 


f^  <0 

n,Tn 


(3.2) 


’n-l/2,m 

k 

^n+l/2,m 


f’*^  >0 

n,m  — 


f*^  <0 

n,in 


For  the  linear  case  [(u.)  =  (po)  =1.0,  K  ,K  ,u,v, 

'^l  m  *^2  n  02 

and  d  constant],  the  constituent  intermediate  time  level  in  the  multi- 


operational  approach,  may  be  eliminated.  The  difference  equation 


obtained  is  consistent  with  the  linear  transport  equation  and  equals  the 
above  difference  equations  plus  some  higher  order  in  time  factorization 
terms.  For  the  nonlinear  case  considered  here,  it  is  not  possible  to 
eliminate  the  constituent  intermediate  time  level.  Therefore  the  exact 
form  of  the  factorization  terms  may  not  be  determined.  However,  their 
numerical  effect  may  be  assessed. 

This  schtoe  is  similar  to  the  standard  ADI  technique  except  that 
upwind  differencing  is  employed  for  the  advective  terms.  The  necessary 
modifications  for  the  X-Sweep  are  shown  in  Table  VII  while  those  em¬ 
ployed  for  the  Y-Sweep  are  given  in  Table  VIII. 

Thus  the  FTUS  scheme  may  be  obtained  from  the  FTCS  scheme  program¬ 
ming  with  only  modest  programming  modification.  In  CRAY-I  FORTRAN  three 
vector  functions  may  be  employed  to  define  the  FTUS  modifications  as 
follows : 


(3.3) 


AMAXI  (x^,  X2)  = 

Xi  >  X2 

(3.4) 

Xz 

Xi  <  X2 

AMINI  (Xj^,  X2)  =  x^ 

Xi  <  Xz 

(3.5) 

Xz 

Xi  >  Xz 

These  functions  eliminate  the  need  for  IF  type  statements. 


CVMGP  (Xj^,  Xz,  x^)  =  Xj^ 

X,  ^  0 

Xz 

*3  <  0 

Table  VIII.  Y-Sweep  Modifications  FTUS 


Equation 

2.11 


2.12 


2.13 


2.13 


2.1A 


2.1A 


FTCS 


FTUS 


k+1 

^n-l/2,m 

2 


max 


k+1 

^n-l/2,m 


V 


k+1 

n+l/2,m 


2 


min 


k+1 

^n+l/2,m 


k+1/ 2* 
n,m-l/2 


1 

-^+1/2* 

n,m+l/2 


1 

-^+1/2* 
n,in+l/ 2 


k+1/ 2* 
s 

n,m 

k+1/ 2* 
^n,m+l 


k+1/ 2* 
'^n,m+l/2  — 


0 


k+1/ 2 
%,m+l/2 


k+1/2* 

n,m-l/2 


-^+1/2*  k+1/2* 

n,m-l/2  ^n,m-l 


-jk+1/2*  k+1/2* 

^n,m-l/2  ^n,m 


k+1/2* 
^n,m-l/2  — 


0 


k+1/2* 

'^n,m-l/2 


<  0 


I 


4 .  The  Spread  Time  Derivative  STCS  Scheme 

Using  previous  notation,  the  approximations  for  the  X-Sweep  may 
now  be  written  as  follows 


f6^(ds)  +  ^  ^  (ds)*^  ^  (l.-) 

1  1  in  ^  - 


At  6 


2Aai(Ui)^ 


a,  6  ,  .k+1/2* 

^+1/2*  „k+l/2* 

a  K  - ; - .  ■  .  —  — 

“l 


•^N-1 


(4 


AC  6 

_ !2. 

2Aa^(M2)„ 


«  (sS 

_ 

(''2)n. 


0 


at  (n,m) 


If  we  place  all  terms  at  time  level  k+1/2*  on  the  left-hand  side  of 
the  equation  and  expand  we  obtain  (2.3)  if 


2 

3 


(ds) 


k+1/2* 

n,m 


+ 


(ds) 


k+1/2* 

n.nri-l 


+  (ds) 

6 


k+1/2* 

n,m-l 


(ds) 


k+1/2* 

n,m 


(4.2) 


and  (2.4)  for 


2 

3 


n.  + 

n  ,m 


(ds) 


k 

n+l,m 


+  (ds) 

6 


k 

n-l,m 


(ds) 


n,m 


(4.3) 


All  other  terms  in  (2.3)  and  (2.4)  remain  the  same.  Equation  (4.2) 
necessitates  the  following  modifications  to  (2. 5)- (2. 7). 


k+1/2* 

a  =  a  + 

n,nH-l  n,m+l  6 


(A. 5) 


,k+l/2* 


n,m 

a  =  a - ^ 

n,m  n,ni  3 


(A. 6) 


The  approximations  for  the  Y-Sweep  are  as  follows 


.  ..k+l  .k+1  At  6  /a„  a-\ 

,k+l/2*,  “y'''  *  ,  “2  (-2.-2-)  k+l  k+l 

3 - 5 -  2to.(V:)  2  '' 

z  z  n  *■ 

oi„  6  (s^"*^^')  At  6  /a.  a  \ 

“2  ^+1  „k+l  “2  .  “1  l-A— )  k+1/2*  k+1/2* 

-  d  K  - 7 - r -  +  -tt: - 7 — ^  \  ds  /  u 

,)  “2  ^^2^n  2Aa^(pp^L  J 

z  n  ^  j 


A4.  i  X  /  k+1/2*. 

At  6  a.,  6  (s  ) 

-ik+1/2.  ^k+1/2. 


If  we  place  all  terms  at  time  level  k+l  on  the  left-hand  side  of  (A. 7) 
and  expand  we  obtain  relations  similar  to  (2. ll)-(2. lA) . 

In  fact  for  (2. ll)-(2. 13) 


n-1  ,m 


a  -  -  —  a  - 

n-l,ni  n*l,m 


(A. 8) 


n+1 .  m 


In  Equation  (2.14),  we  employ  (4.2) 


(ds) 


k+1/2*  ^ 
n,m 


2(ds) 


k+1/2* 

n,m 


3 


+ 


(ds) 


k+1/2* 

n,m+l 


+  (ds) 
6 


k+1/2* 

n,m-l 


(4.11) 


We  therefore  note  that  the  spread  time  derivative  scheme  may  be 
obtained  from  the  standard  scheme  with  only  minor  modifications. 

5 .  Flux-Corrected  Transport 

As  in  the  linear  case,  both  higher  and  lower  order  in  space 
schemes  are  employed.  For  the  nonlinear  case,  the  following  flux  format 
is  needed. 

(5.1) 

^n,m-l/2 


d^^^s^  =  d^  s^  -  rAa^(li^)  Aa„(v„)  1  /o 

n,m  n,m  n,m  n,m  Lllm22nJ  V  n+1/2 


»in 


^n-l/2,m  ^  ^n,m+l/2 


where  t  =  kAt  ,  x  =  2  ^  “X  ^’^2^i^“2 

E  concentration  at  location  (n,m)  at  time  level  k 

n,m 

Aai(Pi)  E  X  space  step  at  m 
JL  1  m 


Aao(Po)  =  y  space  step  at  n 

I  E  general  index  at  time  level  k+1  ,  which  we  set  to 
H  or  L  for  the  higher  or  lower  scheme, 
respectively 


^^4-1  /9  trvH  /9  ^  fluxes  through  the  appropriate  cell  faces  of  cell 
n+l/2,ii+l/2  (n,ni).  Form  dependent  upon  the  finite  difference 

formulation. 


We  observe  from  (5.1)  that  the  difference  between  the  higher  and  lower 
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order  scheme  at  (n,m)  may  be  written  as  follows: 


(s^  -  )=  -fAa^Cy  )  Aa,(y,)  [Yf^ 

\  n,m  n,m/  L  1  1  m  2  2  n  n,nu  LV  n 


'  _  p 

n+l/2,m  n+1/2, 


"  (^n-l/2,m  ^n-l/2,m)  (^n,nri-l/2  ^n,m+-l/2) 


\  n, 


—  v‘~‘ 

m-1/2  n,m-l/2> 


Note  this  difference  may  be  expressed  as  an  array  of  fluxes  between 

adjacent  grid  points.  We  next  develop  the  flux  expressions  for  the 
H  L 

higher  (F  )  and  lower  (F  )  order  schemes.  In  order  to  aid  in  notation, 
we  make  the  following  definition  for  an  arbitrary  variable,  F  . 


A+1/2 


+  F^  NA. 

\  n,m  n,m// 


(5.3) 


For  the  higher  order  scheme  we  employ  the  FTCS  scheme  written  in 

(2.1)  in  which  the  factorization  terms  developed  in  the  multioperational 
method  are  not  shown.  Equation  (2.1)  may  be  written  in  the  form  of 

(5.1) ,  where  the  total  fluxes  are  presented  as  the  sum  of  advective 
and  diffusive  fluxes  as  given  in  Part  II  (7.4)  with  I  -  n  . 

From  Equation  (2.1)  one  then  obtains  for  the  advective  fluxes: 


"A  k+1/2  ,  /  S“  + 

^n+l/2,m  '^n+l/2,m^*^^^l^m^®l  \  2 


S«  + 


,k+l/2 


n+l,m 
n+1 , m  — 


I 

n,m''*“  J/ 


(5.4) 
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“a  k+1/2  .  . 


s”  +  S^  \ 

,  2  )  ^  n.rn+l 
v  /n,in+l  — 


+  (s!L^) 

\  2  /  n.m 

'  'n,ni 


The  diffusive  fluxes  are  then  given  by  the  following  relations. 

(K  E  K  ,  K  5  K  ) 

'  X  ’  y  ct^ 


H 


■0  ^  k+1/2 


n+l/2,m  -  y  . 1 /o  2 

—  ■^n+l/2,m 

_  -  (S“  +  S"^) 


n,m 


,,  1  d. 

n+l,iy  \  n+l,m  n.m  / 


^“2^^2^n+l/2 


(5.5) 


(5.6) 


“o  k+1/2  ^^^^2^n^°‘2 

n.“tl/2  “  -\.:„.i/2  2 


[(s“ 


+  S^)  -  (S^  +  S*^) 

n,m 


^,1  *  d''«^2) 

n,n+lj  \  n,m+l  n.m  / 


^“l^^l^i.H-1/2 


(5.7) 


For  the  lower  order  scheme,  the  FTUS  scheme  written  in  (3.1)  is 
employed.  Factorization  terms  generated  by  the  multioperational  method 
are  not  considered.  Equation  (3.1)  is  written  in  the  form  of  (5.1). 

The  total  fluxes  are  presented  as  the  sum  of  advective  and  diffusive 
fluxes . 

From  Equation  (3.1)  one  obtains  the  following  set  of  advective 
fluxes . 
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rn.l/2.n.^  ^ 

L.  I  n.m 


n+l/2,m 


n-l/2,m 


n,nH-l/2 


n,ni-l/2 


k+l/2  ^  . 

''ti+l/2,m  ° 

^k+1/2 

11+1/2,1 

k+1/2 

1  /o  ="0 

n-l/2,m  “ 

k+1/2 

n-l/2.n 

k+1/2 

Vl/2.m  "  ° 

k+1/2 

n-1/2,1!) 

ls» 

|v 

o 

k+1/2 

“n.i»+l/2‘ 

n,iii+l/2 

k+1/2 

“n,m+l/2‘ 

k+1/2  ^  - 

‘n,m-l/2  — 

k+1/2 
'"n, 01-1/2' 

k+1/2  ,  n 

'n,iii-l/2 

k+1/2 

'*n,m-l/2‘ 

,>^+1/2  ^  +  .k+1/2 

n+l/2.m*‘<''l>m*“l\  2  ,  Vl.m 


The  diffusive  fluxes  are  obtained  from  relations  (5.6)  and  (5.7)  with 
H  replaced  by  L. 

The  Zalesak  flux-correction  procedure  as  reported  In  [3]  for  the 
linear  case  proceeds  analogously  as  follows: 

First,  the  antl-dlf fuslve  fluxes  are  computed. 

^0  ^0 

^n+l/2,m  ^n+l/2,m  ^n+l/2,m  ^  ^n+l/2,m  ^n+l/2,m 


Ha  L.  L. 

-  W  ^  +F^  -F^  (5  13 

^n,nH-l/2  n,nri-l/2  **i,nH-l/2  n,nH-l/2  n,in+l/2  '  ■' 


In  computing  the  difference  between  the  diffusive  fluxes  (third  and 

y 

fourth  terms  in  the  above  expressions) ,  note  that  the  terms  with 
may  be  completely  eliminated. 

The  above  anti-diffusive  fluxes  are  screened  in  the  following 
manner . 


A  /o  =0 
n+l/2,m 


( 

L 

L  \ 

if 

n+l/2,m\ 

.  n+l,m 

-  ^n.m)  • 

and  either 

\+l/2,m( 

^  n+2 ,  m 

- 

n+l,in 

/ 

oL 

_L  \ 

or 

n+l/2,mV 

S 

^n-l,m/ 

(5.14) 


^n,iiri-l/2  ^ 


n,nri-l/2 


and  either  A 


h,nH-l/2 


n,m+l/2 


(s^  )  <  0 

\  n,nH-l  n,m/ 

(s^  -  l)  ^  0 

V  n,m  n,m-l/ 


Next  the  maximum  and  minimum  cell  values  are  determined. 


S  =  max 
n,m 


-max 

S  =  max 
n,m 


(s^  ,s^  )  =  minfs*^  ,S^  ) 

\  n,m  n,m/  n,m  \  n,m  n,m/ 

naxfs^  ,S®  ,S®  ,  ,S®  S®  ) 

\  n-l,m’  n,m’  n+l,m’  n.m-1’  n.m+1/ 


(5.16) 


(5.17a) 


S  =  min 


is;,  .  ,S-  ,S"  T,S"  )  (5.17b  bis) 

\  n-l,m  n,m  n+l,m  n,m-l  n,m+l/  ' 


The  author  has  also  employed  only  quantities  at  time  level  k  ,  obtain¬ 
ing  the  following  alternative  relations: 


max 

S  =  max 
n,m 


=  min(s^  ,S^  ,s\  ,S^  (5.18b) 

n,m  \  n-l,m  n,m  n+l,m  n,m-l  n,m+l/ 


Next  the  sum  of  all  anti-diffusive  fluxes  into  cell  (n,m),  P  , 

’  n,m 


is  determined. 


K,m  =  "'^^<°*Vl/2,m>  -  "‘^"^°’Vl/2.m> 


(5.19) 


+  max(0,A^^^_^^2>  " 


The  maximum  allowable  mass  into  the  cell,  Q  ,  is  then  computed. 

n,m 


+  ^  Amax  _  gL  'jT  )  Aa  (5.20a) 

n,m  \  n,m  n,m/L  1  m  1  2  n  2  n,mj 


Note  S  is  as  given  by  (5.17a).  The  author  has  employed  two  alter- 
n  ,m 

native  formulations. 


+  ^  /'s’"®’'  -  Aa.  (p,)  Aa„d'''*’^l  (5.20b) 

^n.m  V  n.m  n.m/L  1  m  1  *^2  n  2  n.mj 


where  S  is  now  given  by  (5.18a). 
n  ,ni 

The  second  formulation  considered  is  given  by: 


+  ^  Lmax  _  gn  )  Aa^  (p_)  Aa„d’^'''^'|  (5.20c  bis) 

^n,m  \  n,in  n,m/L  1  m  1  '^Z'n  2  n,mj 


where  is  given  by  (5.18a). 

n  yin 

Similarly,  the  sum  of  all  anti-diffusive  fluxes  out  of  cell  (n,m), 

P  ,  is  determined. 
n,m 


P  =  max(0,A  )  -  min(0,A  , ) 

n,m  ’  n+l/2,m  ’  n-l/2,m 


(5.21) 


The  maximum  allowable  mass  to  leave  the  cell,  Q  ,  is  then  computed. 

n,m 


Q  =  (s^  -  s"'^’^^r()i,)  Aa,  (vi„)  Aa„d^''^^1  (5.22a) 

^n,m  V  n,m  n,m/L  1  m  1  2'n  2  n,mj 


Note  S  is  as  given  by  (5.17b). 
n  ,m 

As  in  the  case  of  ,  the  author  has  employed  two  correspond- 

n  y  m 

ing  alternatives.  Under  the  first  alternative,  (5.22a)  is  employed 

with  now  given  by  (5.18b).  Under  the  second  alternative, 

n,m  ® 


Q"  =  (s"  -  s'"^")r(p  )  Aa/y„)  Aa„d'^'^^  1  (5.22b) 

n,m  V  n,m  n,m/L  1  m  l'’^2  n  2  n,mj 


with  S  given  by  (5.18b). 


The  following  ratios  are  next  computed  for  use  In  determining 


■ 

K' 


I 


I 


F.-- 


r  , 


E;;: 


the  limiting  coefficients, 


n,m 


mln^l.Q^  /p^  ^ 

\  n,m/  n,m/ 


P'*'  >0 

n,m 


P"^  =0 

n,m 


(5.23) 


n,m 


mln(l,Q"  /p~  ) 

\  ’^n,m/  n,m/ 


The  limiting  coefficients  are  then  given  by 


P  >  0i 
n,m 


P  =0 
n,m 


(5.24) 


1 

[mln(R^,,  ,R 

1  \n+l,m’  n,m/ 

)  ^n+l/2,m  ^  V 

^n+l/2,m  °  1 

^  min^R^  ,R  ^ 

^  \  n,m  n+l,m; 

1  \+l/2,m‘^°) 

1 

i  min/R'*'  ,R~  '' 

)  \  n.mfl’  n,my 

1  ^n,nifl/2  - 

^n,m+l/2  j 

*  min^R^  ,R  ^ 

^  \  n,m  n.nrt-l/ 

'  \,m+l/2  ^ 

(5.25) 


The  antl-dlf fuslve  fluxes  in  (5.12)  and  (5.13)  are  limited  by 
multiplying  by  the  limiting  coefficients  and  the  solution  is  advanced 
to  the  next  time  level. 
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gk+1  ^  gL  _  (  )  Aa,(y„)  (C  ,  ^  A  . ,  „  ^ 

n,in  n,in  [_  1  1  m  2  2  n  n,mj  n+l/2,in  n+l/2,m 


(5.26) 


-  C 


n-1/2 ,m^n-l/2 ,in  ^  ^n,nri-l/2^n,nH-l/2  *^n,m-l/2^n,m-l/2^ 


k+l  L 

We  observe  that  for  C  . ^  =  C  , , =  0  ,  S  =  S  and  for 

n+l/2,m  n,n^l/2  n,m  n,m 

r  =  r  =  1  f>  c^+l  _  oH 

n+l/2,in  n,n^l/2  ’  ’  n,m  ^n.m 

6.  Additional  Flux-Corrected  Transport  Limiters 

min^  / 

In  conjunction  with  (5.18a),  one  could  insist  S™  =  maxlO.O, 

n,in  \ 

ginin  )  ^fhere  on  the  right  hand  side  is  as  determined  in  (5.18b), 

n,m/  n,m 

Equation  (5.22a)  would  be  employed  for  Q  ,  where  would  be 

n,m  n,m 

replaced  by  . 

n  ,m 


As  another  alternative,  one  could  consider  the  following  relations 

„  „max  ,  f.niin 
for  S  and  S 

n,m  n,m 


„max 

S  =  max 
n,m 


(^n-l,m’^n,m’^n+l,m’\,m-l’^n,mfl’^n,m)  (5.27a) 


k+l\ 


s'"""  =  min(s'^  ,  ,8^^  ,S^  ,  ,S^  ,5^-^  } 

n,m  \  n-l,m  n,m  n+l,m  n,m-l  n,m+l  n,m/ 


(5.27b) 


These  relations  could  be  employed  with  (5.20a)  and  (5.22a). 

As  an  additional  alternative,  one  could  employ  the  second  alter¬ 
nate  limiter  of  the  previous  section;  e.g.,  (5.18a),  (5.18b),  (5.20c) 
and  (5.22b),  on  the  first  time  step  and  the  original  Zalesak  limiter  on 
subsequent  time  steps. 
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Clearly,  many  different  forms  of  the  limiter  are  possible.  We  have 
presented  these  additional  limiters  to  outline  the  direction  of  possible 
future  research. 


PART  IV:  ADDITIONAL  NUMERICAL  CONSIDERATIONS 


In  order  to  develop  a  numerical  model  it  is  necessary  to  develop 
the  necessary  approximations  near  both  open  and  closed  boundaries.  This 
then  leads  to  a  discussion  of  the  tridiagonal  matrix  solution  scheme 
necessary  for  each  sweep. 

The  development  of  flow  field  influence  on  the  effective  disper¬ 
sion  coefficients  is  presented  to  close  the  numerical  model.  This 
closure  is  by  no  means  perfect;  however,  it  is  sufficiently  general  to 
permit  model  calibration.  Additional  approaches  may  be  necessary  to 
Incorporate  wind  effects.  The  approach  followed  here  allows  a  determi¬ 
nation  of  the  anticipated  range  of  physical  dispersion.  In  general, 
the  numerical  scheme  will  also  produce  dispersion.  The  model  should  be 
calibrated  by  adjusting  the  dispersion  coefficients  to  values  within  the 
acceptable  physical  range. 

1.  Approximations  Near  Solid  Boundaries 

In  the  hydrodynamic  equations,  the  convective  acceleration  terms 
and  the  eddy  viscosity  terms  must  be  modified  in  the  vicinity  of  the 
boundaries.  This  is  due  to  the  fact,  that  if  the  standard  differencing 
formulae  at  the  boundary  are  used  points  are  referenced  outside  the  grid. 
No  modifications  are  necessary  for  the  continuity  equation.  Since  the 
transport  equation  is  nothing  more  than  a  constituent  continuity  equation 
we  would  anticipate  no  need  to  modify  the  formula.  This  is  indeed  the 
case,  for  the  difference  formulae  for  continuity  are  cell  centered; 
namely,  fluxes  are  evaluated  at  each  cell  face.  The  fluxes  are  merely 
set  to  zero  in  the  standard  formulae  for  no  flow  conditions  across  the 


:■■■,  % 


appropriate  faces.  Let  us  investigate  this  procedure  in  turn  for  each 
difference  scheme.  In  the  X  or  sweep  each  column  in  the  grid  is 

swept  from  top  to  bottom  stSi-ting  and  ending  at  a  boundary.  In  the  Y 
or  02  sweep  each  row  is  swept  from  left  to  right  again  starting  and 
ending  with  a  boundary. 

Let  us  consider  Equation  (2.8)  of  Part  III,  which  we  rewrite  below 
as  the  standard  equation  for  the  X-a^  sweep  . 


k+1/2*  ^  k+1/2*  ^  k+1/2*  „  ,,  ,, 

a  s  1+a  s  =B  (1.1) 

n,m-l  n,m-l  n,m  n,m  n,m+l  n,m+l  m 


TOP  BOUNDARY 

1 .  FTCS  Scheme 

In  (2.5)  of  Part  III  =  Ky^/2*  ^  q  ^ 

n,m-l/2 

no  mass  transfer  through  the  solid  boundary.  Therefore  a  ,  =  0  . 

n,m-l  ’ 

and  one  obtains 


k+1/2*  ^  k+1/2*  „ 

a  s  +a  _.iS  ,.=B 

n,m  n,m  n,m+l  n,m+l  m 


(1.2) 


2 .  FTUS  Scheme 

Exactly  the  same  conditions  hold  for  this  scheme  and  (1.2)  is 
again  obtained. 

3.  STCS  Scheme 

In  this  scheme  referring  to  (4.4)  of  Part  III  -a  ,  = 

n,m-l 

-a  -  (d^'^^y*)/6  and  a  =  .  For  a  true  land 

n,m-l  \  n,m-l  //  n,m-l  \  n,m-l  // 

(  k+l/2*\  / 

^n  m  1  ~  ^  again  (1.2)  is  obtained.  For  the  case  of 
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' '1  ^ .  'I  B’  ■  .'.'■'V./vri  31  ■. ■. •  ..Ti a  .- 


a  flux  restriction  only  at  a  barrier,  (C'D/  6  ii^  0  and  for  this 
scheme  a  barrier  may  not  form  the  end  of  a  computational  segment. 
BOTTOM  BOUNDARY 

All  statements  previously  made  regarding  the  top  boundary  hold 
directly  if  (n,m-l)  is  replaced  by  (n,m+l).  Equation  (1.1)  now 
becomes 


a  ,  +  a  -  B 

n,m-l  n,m— 1  n,m  n,m  m 


(1.3) 


Let  us  consider  (2.10)  of  Part  III,  which  we  rewrite  below  as  the 
standard  equation  for  the  Y-a2  sweep  . 


k+1  ^  k+1  ^  k+1 

a,  s,  +a  s  +a,,  s.,  =B 

n-l,m  n-l,m  n,m  n,m  n+l,m  n+l,m  n 


(1.4) 


LEFT  BOUNDARY 


1.  FTCS  Scheme 


In  (2.11)  of  Part  III 

n“l/2»m  y  I /o  _ 
n— 1/ ^ ,m 


*  0  ;  i . e . ,  there  is 


no  mass  transfer  through  the  solid  boundary.  Therefore  ^  ®  » 


and  one  obtains 


k+1  ^  k+1 

a  s  +a.,  s,i  =B 

n,m  n,m  n+l,m  n+l,m  n 


(1.5) 


FTUS  Scheme 


a  ,  =0  and  (1.5)  is  again  obtained. 

n-i,m 


STCS  Scheme 


1  ,  =  1/6  and  (1.5)  is  obtained  for  a  land  boundary. 

n-l,m  \n-l,m// 
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RIGHT  BOUNDARY 


All  statements  previously  made  regarding  the  left  boundary  hold 
directly  if  (n-l,m)  is  replaced  by  (n+l,m).  Equation  (1.5)  now  becomes 


k+1  ^  k+1 

L  S  I  ci  S 

n-l,m  n— l,m  n,m  n,m 


B 

n 


(1.6) 


2 .  Approximations  Near  Open  Boundaries 

In  the  hydrodynamic  equations,  convective  and  eddy  viscosity  terms 
in  both  motion  equations  must  be  modified  in  the  vicinity  of  these  type 
boundaries.  No  modifications  are  necessary  for  the  continuity  equation 
nor  for  the  transport  equation.  However  fluxes  must  be  specified  across 
the  appropriate  cell  faces.  Let  us  investigate  this  procedure  for  each 
sweep  in  each  scheme. 


TOP  (-)  AND  BOTTOM  (+)  BOUNDARIES 


1.  FTCS  and  FTUS  Schemes 


k+1/2*  k+1/2* 

u  1 and  K  must  be  specified.  In  (2,8)  of  Part  III 


^k+1/2*  .  1  u  • 

S  ,  must  also  be  given. 
n,m-l 


2.  STCS  Scheme 
k+1/2* 


k+1/2*  k+1/2*  -c-  j  T  /^o  a^ 

, ,  K  and  d  ,  must  be  specified.  In  (2.8) 

n,m7l/2  X  ,  n.mTl  ^ 

+  n,n^l/2  + 


of  Part  III  S 


k+1/2* 


n,ni-l  niust  also  be  given. 

As  a  result,  Equations  (1.2)  and  (1.3)  are  again  obtained. 
LEFT  (-)  AND  RIGHT  (+)  BOUNDARIES 


1.  FTCS  and  FTUS  Schemes 


k+1 


n-l/2,m 


k+1 

and  K  must  be  specified. 

^n-l/2,m 


In  (2.10)  of 


Part  III  s  ,  must  also  be  given. 
n-l,m 

2.  STCS  Scheme 

If4-1  k4-1  k+1 

V  1  /-i  >  K  »  snd  d  ,  must  be  specified.  In  (2.10) 

"+1/2."  >'„-l/2.n,  V-" 

k+l 

of  Part  III  s  1  must  also  be  given.  As  a  result.  Equations  (1.5) 
n-l,m 

and  (1.6)  are  again  obtained. 

The  specification  of  s  deserves  further  attention.  In  the 
general  case,  s  must  be  given  as  a  function  of  time  over  the  period 
of  concern.  However,  in  estuarine  type  (oscillating)  flows,  it  is  pos¬ 
sible  to  compute  s  based  upon  the  values  at  interior  points  during 
ebb  tide.  The  following  procedure  due  to  Leendertse  [1]  is  presented 
subsequently  with  reference  to  Figure  3.  Diffusion  is  allowed  only  on 
one  face  of  the  boundary  cell. 

TOP  (-)  AND  BOTTOM  (+)  BOUNDARIES 


k+1/2*  ^  k  -  k+1/2*  \  ^ 

®n,m-l  ®n,m-l  “n,m-l/2  ^ 


/  k  _  k 
(®n,m  ®n,m-l 


-  ^k+1/2*  I 

''n,m-l/2  ^ 


:y  At 
2  2 


(2.1) 


LEFT  (-)  AND  RIGHT  (+)  BOUNDARIES 


/k+1/2*  _  ^k+l/2*\ 
I  n,m  n-l,m  j 


gk+1  ^  ^k+1/2*  -  ^k+1  \  *  _ +  • 

n-l,m  n-l,m  n-l/2,m  ^“2^‘^2^n-l/2 


/k+1/2*  _  ^k+l/2*\ 
I  n,m  07!, m  / 


(2.2) 


Consider  at  an  arbitrary  boundary  location  b  ,  ebb  flow  conditions 
to  endure  over  N  time  steps.  Compute  at  each  time  step  of  length 
At  the  following  variables. 

Me^  =  Se^  Ve^  Vej  =  At(W)^  (2.3) 

where 

Me^  =  mass  flux  across  boundary  face  during  time  step  i 

Ve^  =  volume  through  boundary  face  during  time  step  i 
o  i  - 

~  concentration  at  the  boundary  face  during  time  step  i 
(Vd)^  =  average  discharge  through  boundary  face  during  time  step  i 
Compute  totals  over  the  ebb  flow  period  as  follows 

N  N 

"ebb  -  2  ""i  ''ebb  ■  E 
i=l  i=l 

where 

^ebb  ”  total  mass  flux  across  boundary  face  during  ebb 
^ebb  -  total  volume  across  boundary  face  during  ebb 
For  the  next  M  time  steps  occurring  during  the  flood  period,  compute 
the  following  quantities. 


^^1  =  '^ebb  M^l  =  ^ebb 


Vf 


i-1  i-1 

1  ■  ''ebb  -  Z  ■  "ebb  -  2  "'‘N-k  "  i  ^ 

k=l  k=l 


(2.5) 


Then  the  boundary  concentration  during  the  flood  period  is  given  by 


(2.6) 


I 


<  = 


(1  -  . 

Vf . 

1 


+  r,  c,  ,  Vf .  >  0 

b  D  1  — 


Vf.  <  0' 
1 


where 


A 


Sfj^  =  boundary  concentration  during  flood  for  time  step  i 
Mf^  =  mass  in  ebb  storage  remaining  prior  to  time  step  i 
Vf^  =  ebb  volume  remaining  prior  to  time  step  i 
r  =  exchange  ratio  for  the  boundary  (0  <  r,  <  1) 
c^  E  ocean  background  concentration 
The  method  above  represents  one  approach  toward  specifying  s  during 
the  flood  period.  Other  approaches  similar  in  concept  are  available. 
Usually  one  must  select  individually  the  most  appropriate  approach  for 
each  area  to  be  modelled. 

3.  Tridiagonal  Matrix  Solution 

Consider  the  following  system  of  equations,  which  is  termed  a 
trldlagonal  system.  We  note  from  our  previous  sections  on  boundary  con¬ 
ditions  that  this  system  of  equations  is  obtained. 


b^vi  +  c^V2 


a3V2  +  b3V3  +  C3V^ 


a.v.  ,  +  b.v.  +  c  .v.^ , 
i  1-1  1  i  1  1+1 


=  d. 


=  d. 


=  d. 


=  d. 


(3.1) 


^N-l’^N-2  ’^N-l'^N-l  '^N-l'^N  ‘^N-1 


N  N-1  N  N  N 
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A  numerically  effective  approach  in  solving  this  system  is  forward 
elimination  and  backward  substitution.  Each  step  is  discussed  in  turn. 
Forward  elimination 


(3.2) 


Substituting  (3.2)  into  the  second  equation  in  (3.1) 


Let  us  define  6^  -  bj^  and  then  we  note  (3.2)  may  be 

written  as  ~  and  (3.3)  may  be  written  as 

(d2  -  ®2^1^  “  ''2''3  ^2 

"2  =  -  "2  =  ^2  -  ^^3 

4  -  ^2  b^  ' 

(3.4) 

""l  ^^2  ~  Vl 

“*2  '  ^2  -  -2  ’'2  ■ 

Let  us  suppose  that  for  equation  i  ,  we  have  the  following  quantities 
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For  the  (i+1)  equation  we  have 


Substituting  the  above  relation  for  v  , 


into  equation  N  we  obtain 


(3.9) 


~  Vn-1 


Y 


N 


where  6  =b  -(ac  This  completes  the  forward  elimination 

N  N  N-1  N-1 

step.  The  are  computed  for  i  =  1,  N  as  summarized  below. 


Si  =  bf  - 


^i^i-1 

^i-1 


(3.10) 


^i  = 


d.  - 
1 


^^i-1 


'i 


=  2,  N 


Backward  substitution 


We  note  v  =  y  and  it  is  therefore  possible  to  then  employ 
N  N 


'i-1 


"i-1  “  "i-1  -  bTT  "i  '  = 


(3.11) 


^i-1 


thereby  computing  employing  the  previously  computed  values 

1  ’  l"’’^!  ’  ^1  forward  elimination  step. 

Let  us  make  the  following  variable  assignments  in  anticipation  of 

coding  efforts. 

O  E  V 


i. 


Thus  we  may  reformulate  our  solution  procedures  as  indicated  in  Table  IX. 


Mitchell  and  Griffiths  [5]  report  that  in  order  for  this  algorithm  to 
be  implemented  on  a  digital  computer  the  following  characteristics  must 
be  possessed  by  the  coefficient  matrix  elements. 


a.  <  0  ,  b.  >  0  ,  c.  <  0  i  =  1,  N 
111 


(3.13) 


b  .  >  -  (a .  +  c . ) 
1  —  1  1 


In  this  case  |P^|  ^1  for  i  =  1,N  and  the  growth  of  errors 
will  be  eliminated.  The  values  of  3  ,  P  ,  and  Q  may  be  output  dur¬ 
ing  preliminary  computations  to  check  the  numerical  formulation. 

4.  Hydrodynamic  Interface 

Within  the  hydrodynamics  program,  a  three  time-level  stabilizing 
correction  scheme  is  employed  to  compute  the  field  variables  n  >  u  , 
and  V  at  each  time  step.  Within  the  two  time-level  transport  scheme, 
u  and  V  are  employed  at  k  ,  k+1/2*  ,  and  k+1  .  Two  approaches 
suggest  themselves  for  interfacing  the  two  schemes. 

Approach  1 

1.  Employ  At  in  the  transport  scheme  equal  to  At„  in  the 

1  n 

hydrodynamic  scheme.  Define  +  u^)/2.  and 

k+1/2*  ,  k+1  ^  k, 

V  =  (v  +  V  )/2.  . 

2.  Perform  one  sweep  in  the  transport  scheme  for  each  sweep  in 


the  hydrodynamic  scheme. 


Approach  2 


1.  Employ  At_  in  the  transport  scheme  as  twice  At  in  the 
1  K 
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Table  IX.  Tridiagonal  Matrix  Solution  Procedure 


Forward  Elimination 


Backward  Substitution 


^^2  - 


^2^ 


> 


®2 

d2  -  a2Qi 


'^N-1  ^N-1  ^N-l'^N 


Vi  =  Qi  - 


■  ®k^k-l 


>  k  =  i,...N 


\  -  ^'^k-l 


K  =  5.93  h  u*  (5.1) 

X 

where 

=  longitudinal  dispersion  coefficient 
h  s  water  depth  (hydraulic  radius) 
u*  =  shear  velocity 

For  open  channel  flow  u*  =  /ghS^  and  from  the  Chezy  relation 

u  =  c/hS  .  As  a  result,  we  obtain 
e 

(5.2) 

where 
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u  r  velocity 
g  E  gravity 
c  =  Chezy  coefficient 
Therefore,  Equation  (5.1)  becomes 

K  =  5.93/g  —  (5.3) 

Taylor  [7]  has  conducted  pipe  flow  experiments  to  determine  the 
longitudinal  dispersion  coefficient.  By  assuming  the  hydraulic  radius 
as  half  the  pipe  radius  in  the  pipe  experiments  and  equal  to  the  water 
depth  in  a  uniform  steady  flow  open  channel,  the  coefficient  in  (5.1) 
and  (5.3)  is  determined  to  be  20.2  rather  than  5.93.  As  a  result,  we 
would  expect  a  general  relationship  of  the  following  form  to  hold. 

=  Cx  e  (5.93,  20.2)  (5.4) 

Wind  and  wave  effects  will  increase  the  effective  dispersion  coeffi¬ 
cients.  The  relationships  are  not  well  known.  However,  Swain  [8]  has 
suggested  the  addition  of  the  following  term  to  (5.4)  to  account  for  wind 
effects 


KW  E  wind  effect  addition 

X 

8  E  ratio  of  sediment  mass  transfer  coefficient  (e  )  to  turbu- 

s 

lent  transfer  coefficient  (e  )  (1  <  6  <  5) 

m  —  — 

u^  E  wind  velocity 
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c  =  Chezy  coefficient 


h  H  water  depth 

K  E  von  Karman's  constant  (.41) 

Elder  [6]  has  reported  a  lateral  dispersion  coefficient  similar 
to  (5.1)  with  5.93  replaced  by  0.23.  For  a  general  model  formulation, 
the  following  relations  will  be  initially  considered. 


.  (u  )  h 

K  =  c  +  cw-^!^+  K' 

X  X  C  c  X 


,  (u  )  h 

K  =  c  —  +  cw  — ^  ^  +  K' 

y  y  c  c  y 


(5.6) 


(5.7) 


where  c  ,  c  e  (5.93,  20.2)  or  0.23  and  are  spatially  variable,  cw 
X  y 

equals  (k»^  6/6)  ,  K'  ,  K'  are  additional  constants. 

X  y 

The  above  relations  represent  a  first  approach  toward  describing 
the  dispersion  mechanisms.  Additional  approaches  may  need  to  be  con¬ 
sidered  to  develop  suitable  results  in  Mississippi  Sound. 
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